In [441]:
import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
from torch.utils.data import Dataset
from torch.utils.data import DataLoader

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from tqdm import tqdm
from time import sleep
import os
from os import listdir
import matplotlib.image as mpimg
from matplotlib.animation import FuncAnimation as FA
import random

pd.set_option('display.max_columns', 50)
pd.set_option('display.max_rows', 50)
In [442]:
# Set Parameters for Training
import matplotlib as mpl
mpl.rcParams['font.family'] = 'serif'  # Use serif font for better readability
mpl.rcParams['font.size'] = 12  # Consistent font size
mpl.rcParams['axes.linewidth'] = 1  # Thicker axes lines
mpl.rcParams['xtick.major.size'] = 5  # Ticks adjustments
mpl.rcParams['ytick.major.size'] = 5
mpl.rcParams['xtick.major.pad'] = 5
mpl.rcParams['ytick.major.pad'] = 5
mpl.rcParams['axes.labelsize'] = 14
mpl.rcParams['axes.titlesize'] = 16
mpl.rcParams['legend.fontsize'] = 12
mpl.rcParams['figure.dpi'] = 300  # High DPI for publication quality
mpl.rcParams['savefig.bbox'] = 'tight'  # Tight bounding box
mpl.rcParams['axes.grid'] = True  # Enable grid by default
mpl.rcParams['grid.alpha'] = 0.5 # Grid alpha
mpl.rcParams['lines.linewidth'] = 1.5 # Line width
mpl.rcParams['scatter.marker'] = 'o' # Marker

Setting seeds for ensuring reproducibility

In [443]:
# import time

# # Generate a random seed based on the current time
# seed = int(time.time())
# print("Selected seed:", seed)

seed = 42
torch.manual_seed(seed)
torch.cuda.manual_seed_all(seed)
np.random.seed(seed)
random.seed(seed)

📌 Loading the data:

There are 100 trajectories of engine degradation in EACH of the datasets (train and test).
For the train data, each engine is run from a certain normal condition till failure.
For the test data, each engine is NOT run until failure (or we have data until a specific point of the engine's state). The amount of time cycles left for this engine to still be 'normal' is what we call RUL, or Remaining Useful Lifetime.

Predicting what the RUL is for the last state of each machine in the test set will be our prediction task.

In [444]:
folder_path = './CMAPSSData/'

listdir(folder_path)

file_name = 'FD003.txt'

df_train = pd.read_csv(folder_path + 'train_' + file_name, header = None, sep = ' ')
df_test = pd.read_csv(folder_path + 'test_'+file_name, header = None, sep = ' ')
rul_test = pd.read_csv(folder_path + 'RUL_'+file_name, header = None)

for df in [df_train, df_test, rul_test]:
    display(df.head())
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
0 1 1 -0.0005 0.0004 100.0 518.67 642.36 1583.23 1396.84 14.62 21.61 553.97 2387.96 9062.17 1.3 47.30 522.31 2388.01 8145.32 8.4246 0.03 391 2388 100.0 39.11 23.3537 NaN NaN
1 1 2 0.0008 -0.0003 100.0 518.67 642.50 1584.69 1396.89 14.62 21.61 554.55 2388.00 9061.78 1.3 47.23 522.42 2388.03 8152.85 8.4403 0.03 392 2388 100.0 38.99 23.4491 NaN NaN
2 1 3 -0.0014 -0.0002 100.0 518.67 642.18 1582.35 1405.61 14.62 21.61 554.43 2388.03 9070.23 1.3 47.22 522.03 2388.00 8150.17 8.3901 0.03 391 2388 100.0 38.85 23.3669 NaN NaN
3 1 4 -0.0020 0.0001 100.0 518.67 642.92 1585.61 1392.27 14.62 21.61 555.21 2388.00 9064.57 1.3 47.24 522.49 2388.08 8146.56 8.3878 0.03 392 2388 100.0 38.96 23.2951 NaN NaN
4 1 5 0.0016 0.0000 100.0 518.67 641.68 1588.63 1397.65 14.62 21.61 554.74 2388.04 9076.14 1.3 47.15 522.58 2388.03 8147.80 8.3869 0.03 392 2388 100.0 39.14 23.4583 NaN NaN
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
0 1 1 -0.0017 -0.0004 100.0 518.67 641.94 1581.93 1396.93 14.62 21.58 554.56 2387.93 9048.65 1.3 47.09 521.89 2387.94 8133.48 8.3760 0.03 391 2388 100.0 39.07 23.4468 NaN NaN
1 1 2 0.0006 -0.0002 100.0 518.67 642.02 1584.86 1398.90 14.62 21.58 554.10 2387.94 9046.53 1.3 47.08 521.85 2388.01 8137.44 8.4062 0.03 391 2388 100.0 39.04 23.4807 NaN NaN
2 1 3 0.0014 -0.0003 100.0 518.67 641.68 1581.78 1391.92 14.62 21.58 554.41 2387.97 9054.92 1.3 47.15 522.10 2387.94 8138.25 8.3553 0.03 391 2388 100.0 39.10 23.4244 NaN NaN
3 1 4 0.0027 0.0001 100.0 518.67 642.20 1584.53 1395.34 14.62 21.59 554.58 2387.94 9055.04 1.3 47.26 522.45 2387.96 8137.07 8.3709 0.03 392 2388 100.0 38.97 23.4782 NaN NaN
4 1 5 -0.0001 0.0001 100.0 518.67 642.46 1589.03 1395.86 14.62 21.58 554.16 2388.01 9048.59 1.3 46.94 521.91 2387.97 8134.20 8.4146 0.03 391 2388 100.0 39.09 23.3950 NaN NaN
0
0 44
1 51
2 27
3 120
4 101

📌 Attaching column names: We have three operational setting columns (os + number), and 21 sensor columns (s + number). We have dropped the last two columns containing NaNs.

In [445]:
col_names = []

col_names.append('unit')
col_names.append('time')

for i in range(1,4):
    col_names.append('os'+str(i))
for i in range(1,22):
    col_names.append('s'+str(i))

df_train = df_train.iloc[:,:-2].copy()
df_train.columns = col_names
display(df_train.head())

df_test = df_test.iloc[:,:-2].copy()
df_test.columns = col_names
display(df_test.head())
unit time os1 os2 os3 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21
0 1 1 -0.0005 0.0004 100.0 518.67 642.36 1583.23 1396.84 14.62 21.61 553.97 2387.96 9062.17 1.3 47.30 522.31 2388.01 8145.32 8.4246 0.03 391 2388 100.0 39.11 23.3537
1 1 2 0.0008 -0.0003 100.0 518.67 642.50 1584.69 1396.89 14.62 21.61 554.55 2388.00 9061.78 1.3 47.23 522.42 2388.03 8152.85 8.4403 0.03 392 2388 100.0 38.99 23.4491
2 1 3 -0.0014 -0.0002 100.0 518.67 642.18 1582.35 1405.61 14.62 21.61 554.43 2388.03 9070.23 1.3 47.22 522.03 2388.00 8150.17 8.3901 0.03 391 2388 100.0 38.85 23.3669
3 1 4 -0.0020 0.0001 100.0 518.67 642.92 1585.61 1392.27 14.62 21.61 555.21 2388.00 9064.57 1.3 47.24 522.49 2388.08 8146.56 8.3878 0.03 392 2388 100.0 38.96 23.2951
4 1 5 0.0016 0.0000 100.0 518.67 641.68 1588.63 1397.65 14.62 21.61 554.74 2388.04 9076.14 1.3 47.15 522.58 2388.03 8147.80 8.3869 0.03 392 2388 100.0 39.14 23.4583
unit time os1 os2 os3 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21
0 1 1 -0.0017 -0.0004 100.0 518.67 641.94 1581.93 1396.93 14.62 21.58 554.56 2387.93 9048.65 1.3 47.09 521.89 2387.94 8133.48 8.3760 0.03 391 2388 100.0 39.07 23.4468
1 1 2 0.0006 -0.0002 100.0 518.67 642.02 1584.86 1398.90 14.62 21.58 554.10 2387.94 9046.53 1.3 47.08 521.85 2388.01 8137.44 8.4062 0.03 391 2388 100.0 39.04 23.4807
2 1 3 0.0014 -0.0003 100.0 518.67 641.68 1581.78 1391.92 14.62 21.58 554.41 2387.97 9054.92 1.3 47.15 522.10 2387.94 8138.25 8.3553 0.03 391 2388 100.0 39.10 23.4244
3 1 4 0.0027 0.0001 100.0 518.67 642.20 1584.53 1395.34 14.62 21.59 554.58 2387.94 9055.04 1.3 47.26 522.45 2387.96 8137.07 8.3709 0.03 392 2388 100.0 38.97 23.4782
4 1 5 -0.0001 0.0001 100.0 518.67 642.46 1589.03 1395.86 14.62 21.58 554.16 2388.01 9048.59 1.3 46.94 521.91 2387.97 8134.20 8.4146 0.03 391 2388 100.0 39.09 23.3950

📌 Attaching RUL(remaining useful lifetime) values to the datasets.

For the train data, the RUL values are not specified, but the document regarding this dataset specifies that for the training data, all the engines were run to failure. Thus, for example, if we had five rows for a specific unit, say unit 7:

unit time cycle
7 1
7 2
7 3
7 4
7 5

Then we know that the last row is when the RUL value becomes 0 (failure), so the RUL for this unit would be attached in this way:

unit time cycle RUL
7 1 4
7 2 3
7 3 2
7 4 1
7 5 0

For the test data, we have the 'solutions' for the test engines in a separate file called rul_test (the name that I used for dataframe). I will use that to attach the RUL values for the test as well.

In [446]:
units_training = max(df_train['unit'])
units_testing = max(df_test['unit'])
print(f'Units in the training dataset: {units_training}')
print(f'Units in the testing dataset: {units_testing}')
df_train = df_train[df_train['unit'] != 260]
if units_training > units_testing:
    df_train = df_train[df_train['unit'] <= units_testing]
elif units_training < units_testing:
    df_test = df_test[df_test['unit']  <= units_training]

assert max(df_train['unit']) == max(df_test['unit'])
no_units = max(df_train['unit'])
print(f'No of units in training and testing dataset after: {no_units}')
Units in the training dataset: 100
Units in the testing dataset: 100
No of units in training and testing dataset after: 100
In [447]:
MAX_RUL = 125
no_units = min(rul_test.shape[0],max(df_train['unit']))
print(f'units :{no_units}')
units :100
In [448]:
rul_list = []
engine_numbers = no_units
for n in np.arange(1,engine_numbers+1):
    
    time_list = np.array(df_train[df_train['unit'] == n]['time'])
    length = len(time_list)
    rul = list(length - time_list)
    rul = [min(MAX_RUL,x) for x in rul]
    rul_list += rul
    
df_train['rul'] = rul_list

rul_list = []

for n in np.arange(1,engine_numbers+1):
    
    time_list = np.array(df_test[df_test['unit'] == n]['time'])
    length = len(time_list)
    rul_val = rul_test.iloc[n-1].item()
    rul = list(length - time_list + rul_val)
    rul = [min(MAX_RUL,x) for x in rul]
    rul_list += rul

df_test['rul'] = rul_list

for df in [df_train, df_test]:
    display(df.head())
unit time os1 os2 os3 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21 rul
0 1 1 -0.0005 0.0004 100.0 518.67 642.36 1583.23 1396.84 14.62 21.61 553.97 2387.96 9062.17 1.3 47.30 522.31 2388.01 8145.32 8.4246 0.03 391 2388 100.0 39.11 23.3537 125
1 1 2 0.0008 -0.0003 100.0 518.67 642.50 1584.69 1396.89 14.62 21.61 554.55 2388.00 9061.78 1.3 47.23 522.42 2388.03 8152.85 8.4403 0.03 392 2388 100.0 38.99 23.4491 125
2 1 3 -0.0014 -0.0002 100.0 518.67 642.18 1582.35 1405.61 14.62 21.61 554.43 2388.03 9070.23 1.3 47.22 522.03 2388.00 8150.17 8.3901 0.03 391 2388 100.0 38.85 23.3669 125
3 1 4 -0.0020 0.0001 100.0 518.67 642.92 1585.61 1392.27 14.62 21.61 555.21 2388.00 9064.57 1.3 47.24 522.49 2388.08 8146.56 8.3878 0.03 392 2388 100.0 38.96 23.2951 125
4 1 5 0.0016 0.0000 100.0 518.67 641.68 1588.63 1397.65 14.62 21.61 554.74 2388.04 9076.14 1.3 47.15 522.58 2388.03 8147.80 8.3869 0.03 392 2388 100.0 39.14 23.4583 125
unit time os1 os2 os3 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21 rul
0 1 1 -0.0017 -0.0004 100.0 518.67 641.94 1581.93 1396.93 14.62 21.58 554.56 2387.93 9048.65 1.3 47.09 521.89 2387.94 8133.48 8.3760 0.03 391 2388 100.0 39.07 23.4468 125
1 1 2 0.0006 -0.0002 100.0 518.67 642.02 1584.86 1398.90 14.62 21.58 554.10 2387.94 9046.53 1.3 47.08 521.85 2388.01 8137.44 8.4062 0.03 391 2388 100.0 39.04 23.4807 125
2 1 3 0.0014 -0.0003 100.0 518.67 641.68 1581.78 1391.92 14.62 21.58 554.41 2387.97 9054.92 1.3 47.15 522.10 2387.94 8138.25 8.3553 0.03 391 2388 100.0 39.10 23.4244 125
3 1 4 0.0027 0.0001 100.0 518.67 642.20 1584.53 1395.34 14.62 21.59 554.58 2387.94 9055.04 1.3 47.26 522.45 2387.96 8137.07 8.3709 0.03 392 2388 100.0 38.97 23.4782 125
4 1 5 -0.0001 0.0001 100.0 518.67 642.46 1589.03 1395.86 14.62 21.58 554.16 2388.01 9048.59 1.3 46.94 521.91 2387.97 8134.20 8.4146 0.03 391 2388 100.0 39.09 23.3950 125
In [449]:
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

# Initialize KMeans and StandardScaler
kmeans = KMeans(n_clusters=6, random_state=42)
scaler = StandardScaler()

# Apply clustering for specific datasets
if file_name in ['FD002.txt', 'FD004.txt']:
    # Ensure required columns are present
    if {'os1', 'os2', 'os3'}.issubset(df_train.columns):
        # Select operating conditions and scale them
        operating_conditions = df_train[['os1', 'os2', 'os3']]
        scaled_conditions = scaler.fit_transform(operating_conditions)

        # Apply k-means clustering and assign operation_mode
        df_train['operation_mode'] = kmeans.fit_predict(scaled_conditions)
    else:
        raise ValueError("Columns 'os1', 'os2', and 'os3' are missing in the dataset!")
else:
    print(f"No clustering applied for dataset {file_name}.")

if file_name in ['FD002.txt', 'FD004.txt']:
    df_test['operation_mode'] = kmeans.predict(scaler.transform(df_test[['os1', 'os2', 'os3']]))
No clustering applied for dataset FD003.txt.
In [450]:
for df in [df_train, df_test]:
    display(df.head())
unit time os1 os2 os3 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21 rul
0 1 1 -0.0005 0.0004 100.0 518.67 642.36 1583.23 1396.84 14.62 21.61 553.97 2387.96 9062.17 1.3 47.30 522.31 2388.01 8145.32 8.4246 0.03 391 2388 100.0 39.11 23.3537 125
1 1 2 0.0008 -0.0003 100.0 518.67 642.50 1584.69 1396.89 14.62 21.61 554.55 2388.00 9061.78 1.3 47.23 522.42 2388.03 8152.85 8.4403 0.03 392 2388 100.0 38.99 23.4491 125
2 1 3 -0.0014 -0.0002 100.0 518.67 642.18 1582.35 1405.61 14.62 21.61 554.43 2388.03 9070.23 1.3 47.22 522.03 2388.00 8150.17 8.3901 0.03 391 2388 100.0 38.85 23.3669 125
3 1 4 -0.0020 0.0001 100.0 518.67 642.92 1585.61 1392.27 14.62 21.61 555.21 2388.00 9064.57 1.3 47.24 522.49 2388.08 8146.56 8.3878 0.03 392 2388 100.0 38.96 23.2951 125
4 1 5 0.0016 0.0000 100.0 518.67 641.68 1588.63 1397.65 14.62 21.61 554.74 2388.04 9076.14 1.3 47.15 522.58 2388.03 8147.80 8.3869 0.03 392 2388 100.0 39.14 23.4583 125
unit time os1 os2 os3 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21 rul
0 1 1 -0.0017 -0.0004 100.0 518.67 641.94 1581.93 1396.93 14.62 21.58 554.56 2387.93 9048.65 1.3 47.09 521.89 2387.94 8133.48 8.3760 0.03 391 2388 100.0 39.07 23.4468 125
1 1 2 0.0006 -0.0002 100.0 518.67 642.02 1584.86 1398.90 14.62 21.58 554.10 2387.94 9046.53 1.3 47.08 521.85 2388.01 8137.44 8.4062 0.03 391 2388 100.0 39.04 23.4807 125
2 1 3 0.0014 -0.0003 100.0 518.67 641.68 1581.78 1391.92 14.62 21.58 554.41 2387.97 9054.92 1.3 47.15 522.10 2387.94 8138.25 8.3553 0.03 391 2388 100.0 39.10 23.4244 125
3 1 4 0.0027 0.0001 100.0 518.67 642.20 1584.53 1395.34 14.62 21.59 554.58 2387.94 9055.04 1.3 47.26 522.45 2387.96 8137.07 8.3709 0.03 392 2388 100.0 38.97 23.4782 125
4 1 5 -0.0001 0.0001 100.0 518.67 642.46 1589.03 1395.86 14.62 21.58 554.16 2388.01 9048.59 1.3 46.94 521.91 2387.97 8134.20 8.4146 0.03 391 2388 100.0 39.09 23.3950 125

Standard normalization (z-score normalization), the data is scaled to have a mean of 0 and a standard deviation of 1.

In [451]:
sensor_colums = ['s'+str(i) for i in range(1,22)]
print(sensor_colums)
['s1', 's2', 's3', 's4', 's5', 's6', 's7', 's8', 's9', 's10', 's11', 's12', 's13', 's14', 's15', 's16', 's17', 's18', 's19', 's20', 's21']
In [452]:
##  Z-score Normalization
if file_name in ['FD001.txt', 'FD003.txt']:
    # Dictionary to store mean and std for columns containing 's'
    mean_std_dict = {}

    # Calculate mean and std for 's' columns in df_train
    for c in sensor_colums:
        mean_std_dict[c + '_mean'] = df_train[c].mean()
        mean_std_dict[c + '_std'] = df_train[c].std()

    # Apply standard normalization to df_train
    for c in sensor_colums:
            df_train[c] = (df_train[c] - mean_std_dict[c + '_mean']) / (mean_std_dict[c + '_std'] + 1e-6)

    # Apply standard normalization to df_test using df_train's mean and std
    for c in sensor_colums:
            df_test[c] = (df_test[c] - mean_std_dict[c + '_mean']) / (mean_std_dict[c + '_std'] + 1e-6)
In [453]:
if file_name in ['FD002.txt', 'FD004.txt']:
    # Function to normalize each group based on its operation_mode
    def z_score_normalize_by_mode(group):
        group = group.copy()  # To avoid modifying the original data
        for c in sensor_colums:
            group[c] = (group[c] - group[c].mean()) / (group[c].std() + 1e-6)
        return group

    # Normalize training data
    df_train_normalized = df_train.groupby('operation_mode', group_keys=False).apply(z_score_normalize_by_mode)

    # Normalize test data using train stats for each operation_mode
    for mode in df_train['operation_mode'].unique():
        train_group = df_train[df_train['operation_mode'] == mode]
        test_group = df_test[df_test['operation_mode'] == mode]
        
        for c in sensor_colums:
            # Ensure the column dtype is compatible
            df_test[c] = df_test[c].astype(float)
                
            mean = train_group[c].mean()
            std = train_group[c].std()
                
            df_test.loc[df_test['operation_mode'] == mode, c] = (
                df_test.loc[df_test['operation_mode'] == mode, c] - mean
            ) / (std+1e-6)

    df_train = df_train_normalized.copy()
In [454]:
if file_name in ['FD002.txt', 'FD004.txt']:
    df_train = df_train.drop('operation_mode', axis = 1)
    df_test = df_test.drop('operation_mode', axis = 1)

os_colums = ['os'+str(i) for i in range(1,4)]
print(os_colums)

mean_std_dict = {}
# Calculate mean and std for 's' columns in df_train
for c in os_colums:
    mean_std_dict[c + '_mean'] = df_train[c].mean()
    mean_std_dict[c + '_std'] = df_train[c].std()

# Apply standard normalization to df_train
for c in os_colums:
    df_train[c] = (df_train[c] - mean_std_dict[c + '_mean']) / (mean_std_dict[c + '_std'] + 1e-6)

# Apply standard normalization to df_test using df_train's mean and std
for c in os_colums:
    df_test[c] = (df_test[c] - mean_std_dict[c + '_mean']) / (mean_std_dict[c + '_std'] + 1e-6)

# Display the first few rows of both datasets
for df in [df_train, df_test]:
    display(df.head())
['os1', 'os2', 'os3']
unit time os1 os2 os3 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21 rul
0 1 1 -0.217019 1.338555 0.0 0.0 -0.187098 -0.712023 -0.780832 -3.552714e-09 0.78153 -0.341487 -0.704766 -0.097136 -0.35337 -0.385592 -0.227589 -0.389847 0.067685 0.469724 -1.734723e-11 -0.889345 0.0 0.0 0.488007 -0.263504 125
1 1 2 0.375359 -1.033981 0.0 0.0 0.080572 -0.497646 -0.775716 -3.552714e-09 0.78153 -0.172752 -0.452059 -0.116655 -0.35337 -0.618867 -0.193798 -0.263362 0.523935 0.729174 -1.734723e-11 -0.321634 0.0 0.0 0.005819 0.375757 125
2 1 3 -0.627127 -0.695047 0.0 0.0 -0.531245 -0.841237 0.116522 -3.552714e-09 0.78153 -0.207663 -0.262529 0.306261 -0.35337 -0.652192 -0.313602 -0.453089 0.361551 -0.100405 -1.734723e-11 -0.889345 0.0 0.0 -0.556733 -0.175053 125
3 1 4 -0.900532 0.321754 0.0 0.0 0.883582 -0.362558 -1.248438 -3.552714e-09 0.78153 0.019257 -0.452059 0.022982 -0.35337 -0.585542 -0.172294 0.052850 0.142818 -0.138414 -1.734723e-11 -0.321634 0.0 0.0 -0.114728 -0.656174 125
4 1 5 0.739899 -0.017180 0.0 0.0 -1.487209 0.080880 -0.697952 -3.552714e-09 0.78153 -0.117477 -0.199352 0.602053 -0.35337 -0.885467 -0.144647 -0.263362 0.217951 -0.153287 -1.734723e-11 -0.321634 0.0 0.0 0.608554 0.437404 125
unit time os1 os2 os3 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21 rul
0 1 1 -0.763830 -1.372915 0.0 0.0 -0.990107 -0.902907 -0.771623 -3.552714e-09 -0.874373 -0.169843 -0.894297 -0.773803 -0.35337 -1.085417 -0.356608 -0.832544 -0.649712 -0.333415 -1.734723e-11 -0.889345 0.0 0.0 0.327278 0.360345 125
1 1 2 0.284224 -0.695047 0.0 0.0 -0.837153 -0.472684 -0.570051 -3.552714e-09 -0.874373 -0.303667 -0.831120 -0.879907 -0.35337 -1.118742 -0.368896 -0.389847 -0.409771 0.165655 -1.734723e-11 -0.889345 0.0 0.0 0.206731 0.587504 125
2 1 3 0.648764 -1.033981 0.0 0.0 -1.487209 -0.924932 -1.284251 -3.552714e-09 -0.874373 -0.213481 -0.641590 -0.459994 -0.35337 -0.885467 -0.292098 -0.832544 -0.360693 -0.675492 -1.734723e-11 -0.889345 0.0 0.0 0.447825 0.210246 125
3 1 4 1.241142 0.321754 0.0 0.0 -0.493006 -0.521139 -0.934313 -3.552714e-09 -0.322405 -0.164024 -0.831120 -0.453988 -0.35337 -0.518892 -0.184582 -0.706059 -0.432190 -0.417694 -1.734723e-11 -0.321634 0.0 0.0 -0.074545 0.570751 125
4 1 5 -0.034749 0.321754 0.0 0.0 0.004095 0.139613 -0.881107 -3.552714e-09 -0.874373 -0.286212 -0.388882 -0.776806 -0.35337 -1.585292 -0.350465 -0.642816 -0.606086 0.304469 -1.734723e-11 -0.889345 0.0 0.0 0.407642 0.013241 125

Below, we can see an example of all the sensor values of a specific engine in the training set (unit 4), as the engine progresses toward failure.

In [455]:
sample = 4
sample_df = df_train[df_train['unit'] == sample].copy()

# Select only sensor columns
n_sensors = len(sensor_colums)
rows = (n_sensors + 2) // 3  # Calculate rows for the grid

# Create subplots
fig, axes = plt.subplots(rows, 3, figsize=(15, 5 * rows))

# Flatten axes for easier iteration
axes = axes.flatten()

# Plot each sensor's data
for j, c in enumerate(sensor_colums):
    axes[j].plot(sample_df['time'], sample_df[c], label='Sensor data')
    axes[j].set_title(c)
    # axes[j].set_xlabel('Time')
    # axes[j].set_ylabel('Value')
    axes[j].legend()

# Hide unused subplots
for ax in axes[n_sensors:]:
    ax.axis('off')

# Adjust layout for better spacing
fig.tight_layout(pad=3.0)
plt.savefig(f'figures\\sensor_of_engine_{sample}_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
No description has been provided for this image

We plot the sensor data obtained from sensor 2 for 5 of the engines to show a general trend.

In [456]:
sensor = 's2'
for sample in range(1,6):
    sample_df = df_train[df_train['unit'] == sample].copy()
    sensordata = sample_df[sensor].to_numpy()
    plt.plot(sensordata, label = "engine "+str(sample))
plt.grid()
plt.legend()
plt.ylabel('Sensor 2')
plt.xlabel('Cycle')
plt.savefig(f'figures\\sensor_{sensor}_of_engines_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
No description has been provided for this image
In [457]:
import seaborn as sns
import matplotlib.pyplot as plt

# Compute correlation matrix
corr_matrix = df_train.corr()

# Plot heatmap
plt.figure(figsize=(20, 14))
sns.heatmap(corr_matrix, annot=True, cmap="coolwarm", fmt=".2f")
plt.title("Feature Correlation Matrix")
plt.savefig(f'figures\\Feature_Correlation_Matrix_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
No description has been provided for this image
In [458]:
from sklearn.ensemble import RandomForestRegressor
import pandas as pd

Xt = df_train.iloc[:,2:-1]
yt = df_train.iloc[:,-1]

rf = RandomForestRegressor()
rf.fit(Xt, yt)

# Get feature importance
importance = pd.Series(rf.feature_importances_, index=Xt.columns)
importance = importance.sort_values(ascending=False)

print("Feature Importances:\n", importance)
Feature Importances:
 s11    0.632116
s9     0.140358
s14    0.031786
s4     0.025796
s12    0.024264
s7     0.020847
s13    0.018195
s8     0.015298
s15    0.013585
s3     0.012977
s21    0.012731
s2     0.012119
s6     0.010617
s20    0.009913
os1    0.009259
os2    0.005749
s17    0.004232
s10    0.000159
s5     0.000000
s1     0.000000
s16    0.000000
s18    0.000000
s19    0.000000
os3    0.000000
dtype: float64
In [459]:
# Combine thresholding and sorting
threshold = 1e-3
low_importance_features = importance[importance < threshold].sort_values()
print("Low importance features sorted:\n", low_importance_features)
Low importance features sorted:
 s5     0.000000
s1     0.000000
s16    0.000000
s18    0.000000
s19    0.000000
os3    0.000000
s10    0.000159
dtype: float64

We can notice here, that there are multiple sensors which are not changing its value. Perhaps, they are not useful features for prediction. Would they have similar behaviors for other engine units as well? Below show that yes (standard deviation is 0, or practically 0).

In [460]:
df_train[low_importance_features.index].describe()
Out[460]:
s5 s1 s16 s18 s19 os3 s10
count 2.472000e+04 24720.0 2.472000e+04 24720.0 24720.0 24720.0 2.472000e+04
mean -3.552714e-09 0.0 -1.734723e-11 0.0 0.0 0.0 -8.342324e-14
std 4.135987e-25 0.0 3.231240e-27 0.0 0.0 0.0 9.997131e-01
min -3.552714e-09 0.0 -1.734723e-11 0.0 0.0 0.0 -3.222112e+00
25% -3.552714e-09 0.0 -1.734723e-11 0.0 0.0 0.0 -3.533705e-01
50% -3.552714e-09 0.0 -1.734723e-11 0.0 0.0 0.0 -3.533705e-01
75% -3.552714e-09 0.0 -1.734723e-11 0.0 0.0 0.0 -3.533705e-01
max -3.552714e-09 0.0 -1.734723e-11 0.0 0.0 0.0 5.384113e+00

📌 Previously mentioned columns are dropped.

In [461]:
#Drop os3, s1, s5, s6, s10, s16, s18, s19 from both train and test
drop_cols1 = low_importance_features.index

df_train = df_train.drop(drop_cols1, axis = 1)
df_test = df_test.drop(drop_cols1, axis = 1)

for df in [df_train, df_test]:
    display(df.head())
unit time os1 os2 s2 s3 s4 s6 s7 s8 s9 s11 s12 s13 s14 s15 s17 s20 s21 rul
0 1 1 -0.217019 1.338555 -0.187098 -0.712023 -0.780832 0.78153 -0.341487 -0.704766 -0.097136 -0.385592 -0.227589 -0.389847 0.067685 0.469724 -0.889345 0.488007 -0.263504 125
1 1 2 0.375359 -1.033981 0.080572 -0.497646 -0.775716 0.78153 -0.172752 -0.452059 -0.116655 -0.618867 -0.193798 -0.263362 0.523935 0.729174 -0.321634 0.005819 0.375757 125
2 1 3 -0.627127 -0.695047 -0.531245 -0.841237 0.116522 0.78153 -0.207663 -0.262529 0.306261 -0.652192 -0.313602 -0.453089 0.361551 -0.100405 -0.889345 -0.556733 -0.175053 125
3 1 4 -0.900532 0.321754 0.883582 -0.362558 -1.248438 0.78153 0.019257 -0.452059 0.022982 -0.585542 -0.172294 0.052850 0.142818 -0.138414 -0.321634 -0.114728 -0.656174 125
4 1 5 0.739899 -0.017180 -1.487209 0.080880 -0.697952 0.78153 -0.117477 -0.199352 0.602053 -0.885467 -0.144647 -0.263362 0.217951 -0.153287 -0.321634 0.608554 0.437404 125
unit time os1 os2 s2 s3 s4 s6 s7 s8 s9 s11 s12 s13 s14 s15 s17 s20 s21 rul
0 1 1 -0.763830 -1.372915 -0.990107 -0.902907 -0.771623 -0.874373 -0.169843 -0.894297 -0.773803 -1.085417 -0.356608 -0.832544 -0.649712 -0.333415 -0.889345 0.327278 0.360345 125
1 1 2 0.284224 -0.695047 -0.837153 -0.472684 -0.570051 -0.874373 -0.303667 -0.831120 -0.879907 -1.118742 -0.368896 -0.389847 -0.409771 0.165655 -0.889345 0.206731 0.587504 125
2 1 3 0.648764 -1.033981 -1.487209 -0.924932 -1.284251 -0.874373 -0.213481 -0.641590 -0.459994 -0.885467 -0.292098 -0.832544 -0.360693 -0.675492 -0.889345 0.447825 0.210246 125
3 1 4 1.241142 0.321754 -0.493006 -0.521139 -0.934313 -0.322405 -0.164024 -0.831120 -0.453988 -0.518892 -0.184582 -0.706059 -0.432190 -0.417694 -0.321634 -0.074545 0.570751 125
4 1 5 -0.034749 0.321754 0.004095 0.139613 -0.881107 -0.874373 -0.286212 -0.388882 -0.776806 -1.585292 -0.350465 -0.642816 -0.606086 0.304469 -0.889345 0.407642 0.013241 125

📌 Splitting Train and Validation Sets: Out of the 100 engines in the training set, I will randomly take out 20 engines for validation.

In [462]:
ratio = 0.8                     #Ratio of training and validation datasets

units = np.arange(1, no_units+1)
no_selected = round(ratio * engine_numbers)
train_units = list(np.random.choice(units,no_selected, replace = False))
val_units = list(set(units) - set(train_units))
print(val_units)


train_data = df_train[df_train['unit'].isin(train_units)].copy()
val_data = df_train[df_train['unit'].isin(val_units)].copy()
[14, 15, 18, 24, 34, 40, 44, 48, 53, 62, 63, 72, 74, 78, 80, 82, 85, 87, 90, 94]

📌 The time series for sensor values were noisy. If the time-series values are $t_1, t_2, t_3, ..., t_n$, then the smoothed values $v_1, v_2, ..., v_n$ with the parameter $\beta$ follow the following formula:

$$v_0 = 0, v_{t} = \frac{\beta v_{t-1} + (1-\beta) x_{t}}{1-\beta^{t}}$$

In [463]:
#Smoothing Function: Exponentially Weighted Averages

def smooth(s, b = 0.98):

    v = np.zeros(len(s)+1) #v_0 is already 0.
    bc = np.zeros(len(s)+1)  

    for i in range(1, len(v)): #v_t = 0.95
        v[i] = (b * v[i-1] + (1-b) * s[i-1]) 
        bc[i] = 1 - b**i

    sm = v[1:] / bc[1:]
    
    return sm

# s = [1,2,3,4,5]
# print(s)
# print(f'After Smoothing: {smooth(s)}')
In [464]:
#Smoothing each time series for each engine in both training and test sets
beta = 0.98
# if file_name=='FD002' or file_name=='FD004':
#     beta = 0.8

# Verify data integrity
assert 'unit' in df_train.columns, "The 'unit' column is missing in df_train"
assert 'unit' in df_test.columns, "The 'unit' column is missing in df_test"

def smooth_series(df, beta, unit_col='unit', sensor_prefix='s'):
    """
    Smooths time series for each engine and sensor column in the dataset.
    """
    for col in df.columns:
        if sensor_prefix in col:
            sm_list = []
            for unit in df[unit_col].unique():
                # Get sensor data for the current unit
                unit_data = df[df[unit_col] == unit]
                s = np.array(unit_data[col].copy())
                sm = list(smooth(s, beta))  # Apply smoothing
                sm_list.extend(sm)  # Append smoothed data
            # Check the length match
            if len(sm_list) != len(df):
                raise ValueError(f"Length mismatch for column {col}: sm_list={len(sm_list)}, df={len(df)}")
            # Add smoothed column to DataFrame
            df[col + '_smoothed'] = sm_list
    return df

# Apply smoothing to training and test sets
df_train = smooth_series(df_train, beta)
df_test = smooth_series(df_test, beta)

Let's take a look at how smoothed values (salmon) look compared to the original series (lightblue), for a particular sensor (unit 10 from train)

In [465]:
df_train.columns
Out[465]:
Index(['unit', 'time', 'os1', 'os2', 's2', 's3', 's4', 's6', 's7', 's8', 's9',
       's11', 's12', 's13', 's14', 's15', 's17', 's20', 's21', 'rul',
       'os1_smoothed', 'os2_smoothed', 's2_smoothed', 's3_smoothed',
       's4_smoothed', 's6_smoothed', 's7_smoothed', 's8_smoothed',
       's9_smoothed', 's11_smoothed', 's12_smoothed', 's13_smoothed',
       's14_smoothed', 's15_smoothed', 's17_smoothed', 's20_smoothed',
       's21_smoothed'],
      dtype='object')
In [466]:
# Set scientific theme
sns.set_theme(style="whitegrid", font="Times New Roman")

# Select unit
sample_unit = 10
sample_df = df_train[df_train['unit'] == sample_unit].copy()

# Extract sensor columns
sensor_cols = [col for col in df_train.columns if col.startswith('s') and '_smoothed' not in col]
n_sensors = len(sensor_cols)

if n_sensors == 0:
    raise ValueError("No valid sensor columns found in the dataset.")

# Define layout for compact visualization
num_cols = 5
num_rows = -(-n_sensors // num_cols)

fig, axes = plt.subplots(num_rows, num_cols, figsize=(9, 2.5 * num_rows), sharex=True, constrained_layout=True)

# Flatten axes
axes = axes.flatten()

# Define colors
colors = {"original": "tab:blue", "smoothed": "#D55E00", "line": "black"}

for j, col in enumerate(sensor_cols):
    axes[j].plot(sample_df['time'], sample_df[col], color=colors["original"], label="Original", lw=0.8)
    
    smoothed_col = f"{col}_smoothed"
    if smoothed_col in sample_df.columns:
        axes[j].plot(sample_df['time'], sample_df[smoothed_col], color=colors["smoothed"], label="Smoothed", lw=0.8)

    # Grid and Labels
    axes[j].set_title(col, fontsize=10)
    axes[j].tick_params(axis='both', labelsize=8)
    axes[j].grid(True, linestyle="--", linewidth=0.5, alpha=0.6)

# Hide unused subplots
for ax in axes[n_sensors:]:
    ax.set_visible(False)

# Global legend for clarity
fig.legend(["Original", "Smoothed"], loc="upper center", ncol=2, bbox_to_anchor=(0.5, -0.02), fontsize=9, frameon=False)

# Save and display
output_dir = "figures"
os.makedirs(output_dir, exist_ok=True)
plt.savefig(os.path.join(output_dir, f"Smoothing_Sensors_{file_name}.png"), dpi=600, bbox_inches="tight", pad_inches=0.05)
plt.show()
No description has been provided for this image
In [467]:
#Remove the original series

for c in df_train.columns:
    if ('s' in c) and ('smoothed' not in c):
        df_train[c] = df_train[c+'_smoothed']
        df_train.drop(c+'_smoothed', axis = 1, inplace = True)
        
for c in df_test.columns:
    if ('s' in c) and ('smoothed' not in c):
        df_test[c] = df_test[c+'_smoothed']
        df_test.drop(c+'_smoothed', axis = 1, inplace = True)
        
for df in [df_train, df_test]:
    display(df.head())
unit time os1 os2 s2 s3 s4 s6 s7 s8 s9 s11 s12 s13 s14 s15 s17 s20 s21 rul
0 1 1 -0.217019 1.338555 -0.187098 -0.712023 -0.780832 0.78153 -0.341487 -0.704766 -0.097136 -0.385592 -0.227589 -0.389847 0.067685 0.469724 -0.889345 0.488007 -0.263504 125
1 1 2 0.082162 0.140305 -0.051911 -0.603752 -0.778248 0.78153 -0.256267 -0.577137 -0.106994 -0.503408 -0.210522 -0.325965 0.298114 0.600759 -0.602622 0.244478 0.059355 125
2 1 3 -0.159060 -0.143790 -0.214927 -0.684518 -0.473946 0.78153 -0.239737 -0.470142 0.033550 -0.554008 -0.245579 -0.369199 0.319689 0.362300 -0.700134 -0.028006 -0.020365 125
3 1 4 -0.350083 -0.023854 0.068078 -0.601573 -0.673476 0.78153 -0.173014 -0.465483 0.030827 -0.562132 -0.226699 -0.260468 0.274122 0.233303 -0.602622 -0.050348 -0.184166 125
4 1 5 -0.123190 -0.022464 -0.255673 -0.459512 -0.678571 0.78153 -0.161453 -0.410085 0.149734 -0.629438 -0.209619 -0.261070 0.262429 0.152830 -0.544131 0.086810 -0.054779 125
unit time os1 os2 s2 s3 s4 s6 s7 s8 s9 s11 s12 s13 s14 s15 s17 s20 s21 rul
0 1 1 -0.763830 -1.372915 -0.990107 -0.902907 -0.771623 -0.874373 -0.169843 -0.894297 -0.773803 -1.085417 -0.356608 -0.832544 -0.649712 -0.333415 -0.889345 0.327278 0.360345 125
1 1 2 -0.234510 -1.030558 -0.912858 -0.685623 -0.669819 -0.874373 -0.237431 -0.862389 -0.827391 -1.102248 -0.362814 -0.608959 -0.528530 -0.081359 -0.889345 0.266396 0.475071 125
2 1 3 0.065883 -1.031722 -1.108189 -0.767009 -0.878781 -0.874373 -0.229286 -0.787298 -0.702443 -1.028523 -0.338765 -0.684998 -0.471450 -0.283418 -0.889345 0.328098 0.385007 125
3 1 4 0.368660 -0.683031 -0.949702 -0.703667 -0.893088 -0.732172 -0.212473 -0.798587 -0.638434 -0.897229 -0.299043 -0.690424 -0.461336 -0.318011 -0.743088 0.224366 0.432860 125
4 1 5 0.284686 -0.473873 -0.751158 -0.528128 -0.890594 -0.761773 -0.227822 -0.713303 -0.667238 -1.040457 -0.309747 -0.680514 -0.491467 -0.188435 -0.773533 0.262517 0.345511 125

📌 When we look at the length of the trajectories for each unit in both the training and test sets, thus 200 engines, we see that for the training set that had the entire trajectories, have minimum length of 128. However, for the test set, we have a minimum value of 31. This means that if we were to predict the final RUL for every test engine unit, we cannot use a window size greater than 31 for training the model. Furthermore, if you take a look at the smoothed series visualization, I put a black vertical line for every graph, at about time step 10. This is because, when we smooth using exponentially weighted averages, seems like in the beginning, it seems a little shaky. Thus, I will be taking that part out of the training process. This will happen for test data as well. Thus, the maximum window size we can take is 31 - 10 = 21. I will go with 20.

In [468]:
print('training set time cycles:')
display(df_train.groupby('unit')['time'].max().describe())
print('test set time cycles:')
display(df_test.groupby('unit')['time'].max().describe())
training set time cycles:
count    100.00000
mean     247.20000
std       86.48384
min      145.00000
25%      189.75000
50%      220.50000
75%      279.75000
max      525.00000
Name: time, dtype: float64
test set time cycles:
count    100.000000
mean     165.960000
std       86.891773
min       38.000000
25%      105.000000
50%      148.000000
75%      208.000000
max      475.000000
Name: time, dtype: float64
In [469]:
n_features = len([c for c in df_train.columns if 's' in c]) + 1 #plus one for time
window = 15
print(f'number of features: {n_features}, window size: {window}')
number of features: 18, window size: 15

📌 Splitting Train and Validation Sets: Out of the 100 engines in the training set, I will randomly take out 20 engines for validation.

In [470]:
train_indices = list(train_data[(train_data['rul'] >= (window - 1)) & (train_data['time'] > 10)].index)
val_indices = list(val_data[(val_data['rul'] >= (window - 1)) & (val_data['time'] > 10)].index)

📌 Normalize the RUL index dividing by maximum value

In [471]:
rul_max = max(df_train['rul'])
df_train['rul'] = df_train['rul'] / rul_max 
df_test['rul'] = df_test['rul'] / rul_max 
rul_max
Out[471]:
125

📌 Z-Normalize the time index

In [472]:
# max_time = max(df_train['time'])
df_train_mean = df_train['time'].mean()
df_train_std = df_train['time'].std()
df_train['time'] = (df_train['time'] - df_train_mean) / df_train_std
df_test['time'] = (df_test['time'] - df_train_mean) / df_train_std
display(df_train)
unit time os1 os2 s2 s3 s4 s6 s7 s8 s9 s11 s12 s13 s14 s15 s17 s20 s21 rul
0 1 -1.396881 -0.217019 1.338555 -0.187098 -0.712023 -0.780832 0.78153 -0.341487 -0.704766 -0.097136 -0.385592 -0.227589 -0.389847 0.067685 0.469724 -0.889345 0.488007 -0.263504 1.000
1 1 -1.386765 0.082162 0.140305 -0.051911 -0.603752 -0.778248 0.78153 -0.256267 -0.577137 -0.106994 -0.503408 -0.210522 -0.325965 0.298114 0.600759 -0.602622 0.244478 0.059355 1.000
2 1 -1.376648 -0.159060 -0.143790 -0.214927 -0.684518 -0.473946 0.78153 -0.239737 -0.470142 0.033550 -0.554008 -0.245579 -0.369199 0.319689 0.362300 -0.700134 -0.028006 -0.020365 1.000
3 1 -1.366531 -0.350083 -0.023854 0.068078 -0.601573 -0.673476 0.78153 -0.173014 -0.465483 0.030827 -0.562132 -0.226699 -0.260468 0.274122 0.233303 -0.602622 -0.050348 -0.184166 1.000
4 1 -1.356414 -0.123190 -0.022464 -0.255673 -0.459512 -0.678571 0.78153 -0.161453 -0.410085 0.149734 -0.629438 -0.209619 -0.261070 0.262429 0.152830 -0.544131 0.086810 -0.054779 1.000
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
24715 100 0.090270 0.142461 -0.166179 1.153957 0.975020 1.250410 0.78153 -0.749628 0.654728 -0.351183 1.280620 -0.719768 0.632469 -0.580150 1.302012 1.086300 -1.285265 -1.267812 0.032
24716 100 0.100387 0.172284 -0.184435 1.163151 0.980149 1.277205 0.78153 -0.753956 0.662013 -0.343280 1.303755 -0.724923 0.646884 -0.567745 1.316765 1.092510 -1.299525 -1.300125 0.024
24717 100 0.110504 0.153569 -0.152430 1.194240 1.010314 1.308290 0.78153 -0.761611 0.671792 -0.336761 1.329882 -0.732999 0.655665 -0.566323 1.332893 1.110514 -1.322759 -1.318544 0.016
24718 100 0.120621 0.128570 -0.121128 1.228660 1.018329 1.328294 0.78153 -0.772824 0.682682 -0.331678 1.346339 -0.737222 0.666908 -0.558203 1.346065 1.116203 -1.352230 -1.335558 0.008
24719 100 0.130737 0.126100 -0.097623 1.250293 1.030726 1.368224 0.78153 -0.782143 0.693334 -0.322193 1.367327 -0.744766 0.675251 -0.555674 1.359017 1.133673 -1.359986 -1.350879 0.000

24720 rows × 20 columns

📌 Prepare Training, Validation and Test Dataloaders.

For training process, I will take batches of 128.
For the validation process, I will take the entire validation set (all the windows avaiable). For the test process, I will take the LAST window of each engine's given trajectory, thus it will have exactly 100 X(size 20 window)s and 100 y(RUL)s.

In [473]:
class data(Dataset):
    
    def __init__(self, list_indices, df_train):
        
        self.indices = list_indices
        self.df_train = df_train
        
    def __len__(self):
        
        return len(self.indices)
    
    def __getitem__(self, idx):
        
        ind = self.indices[idx]
        X_ = self.df_train.iloc[ind : ind + window, :].drop(['unit','rul'], axis = 1).copy().to_numpy()
        y_ = self.df_train.iloc[ind + window - 1]['rul']
        
        return X_, y_
    
train = data(train_indices, df_train)
val = data(val_indices, df_train)

trainloader = DataLoader(train, batch_size = 256, shuffle = True)
valloader = DataLoader(val, batch_size = 128, shuffle = True)

units = np.arange(1, no_units+1)

class test(Dataset):
    
    def __init__(self, units, df_test):
        
        self.units = units
        self.df_test = df_test
        
    def __len__(self):
        
        return len(self.units)
    
    def __getitem__(self, idx):
        
        n = self.units[idx]
        U = self.df_test[self.df_test['unit'] == n].copy()
        X_ = U.reset_index().iloc[-window:,:].drop(['index', 'unit','rul'], axis = 1).copy().to_numpy()
        y_ = U['rul'].min()
        
        return X_, y_
    
test = test(units, df_test)
testloader = DataLoader(test, batch_size = 100)
In [474]:
dataiter = iter(trainloader)
x,y = next(dataiter)
x.shape
Out[474]:
torch.Size([256, 15, 18])

2. Model Building and Training¶

📌 I will be using a Model incoprorating Attention Mechanism with pairwise interactions and positional awareness layers before the final output.

In [475]:
def anderson(f, x0, m=5, lam=1e-4, max_iter=50, tol=1e-2, beta=1.0, verbose=False):
    """Improved Anderson acceleration for fixed-point iteration."""
    # Shape parameters
    batch, channels, dim = x0.shape
    
    # Storage for historical X and F
    X = torch.zeros(batch, m, channels * dim, dtype=x0.dtype, device=x0.device)
    F = torch.zeros(batch, m, channels * dim, dtype=x0.dtype, device=x0.device)

    # Initialize X and F with the first two iterations
    X[:, 0], F[:, 0] = x0.view(batch, -1), f(x0).view(batch, -1)
    X[:, 1], F[:, 1] = F[:, 0], f(F[:, 0].view(batch, channels, dim)).view(batch, -1)

    # Prepare the H matrix and y vector for Anderson acceleration
    H = torch.zeros(batch, m + 1, m + 1, dtype=x0.dtype, device=x0.device)
    H[:, 0, 1:] = H[:, 1:, 0] = 1  # First row and column for normalization
    y = torch.zeros(batch, m + 1, 1, dtype=x0.dtype, device=x0.device)
    y[:, 0, 0] = 1  # y vector for the Anderson step

    res = []  # To store residuals

    for k in range(2, max_iter):
        n = min(k, m)  # Number of iterations to consider

        # Compute residual matrix G
        G = F[:, :n] - X[:, :n]

        # Compute the H matrix (Gramian + regularization)
        GTG = torch.bmm(G, G.transpose(1, 2))
        H[:, 1:n+1, 1:n+1] = GTG + lam * torch.eye(n, dtype=x0.dtype, device=x0.device)[None]

        # Solve for alpha using least squares
        try:
            alpha = torch.linalg.solve(H[:, :n+1, :n+1], y[:, :n+1])[:, 1:n+1, 0]  # (batch x n)
        except RuntimeError as e:
            if verbose:
                print(f"Solver failed at iteration {k}: {e}")
            break

        # Update X and F
        X_update = beta * (alpha[:, None] @ F[:, :n])[:, 0] + (1 - beta) * (alpha[:, None] @ X[:, :n])[:, 0]
        idx = k % m
        X[:, idx] = X_update
        F[:, idx] = f(X[:, idx].view(batch, channels, dim)).view(batch, -1)

        # Compute residuals
        residual = (F[:, idx] - X[:, idx]).norm(dim=1) / (1e-5 + X[:, idx].norm(dim=1).clamp(min=1e-5))
        res.append(residual.mean().item())  # Mean residual across batch

        if verbose:
            print(f"Iteration {k}, residual: {res[-1]:.6f}")

        # Check for convergence
        if res[-1] < tol:
            if verbose:
                print(f"Converged at iteration {k}, residual: {res[-1]:.6f}")
            break

    return X[:, idx].view(batch, channels, dim), (res, k)
In [476]:
import torch.autograd as autograd

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print('Using device:', device)

class DEQFixedPoint(nn.Module):
    def __init__(self, f, solver, **kwargs):
        super().__init__()
        self.f = f
        self.solver = solver
        self.kwargs = kwargs

    def forward(self, x):
        # Compute forward pass and re-engage autograd tape
        b, c, f = x.shape
        with torch.no_grad():  # Ensure no gradients for intermediate calculations
            z, self.stats = self.solver(lambda z: self.f(z, x), torch.zeros((b, 4 * c, f), device=device), **self.kwargs)

        z = self.f(z, x)

        # Set z to require gradients before the hook
        z = z.requires_grad_()  # Ensure z requires gradients

        # Set up Jacobian vector product (without additional forward calls)
        z0 = z.clone().detach().requires_grad_()  # Ensure z0 requires gradients
        f0 = self.f(z0, x)

        def backward_hook(grad):
            # Perform the backward pass using the solver
            g, self.backward_res = self.solver(
                lambda y: autograd.grad(f0, z0, y, retain_graph=True)[0] + grad,
                grad, **self.kwargs)
            return g

        # Register the backward hook
        z.register_hook(backward_hook)
        return z, self.stats
Using device: cuda
In [477]:
import torch.optim.lr_scheduler as lr_scheduler

# class ScaledLeakyReLU(nn.Module):
#     def __init__(self, negative_slope=0.1, scaled_slope=0.9):
#         super(ScaledLeakyReLU, self).__init__()
#         self.negative_slope = negative_slope
#         self.scaled_slope = scaled_slope

#     def forward(self, x):
#         return torch.where(x > 0, self.scaled_slope * x, self.negative_slope * x)

class DualInputAttention(nn.Module):
    def __init__(self, channels, feature_size):
        super(DualInputAttention, self).__init__()
        self.feature_size = feature_size
        self.query_proj = nn.Linear(self.feature_size, self.feature_size)
        self.key_proj = nn.Linear(self.feature_size, self.feature_size)
        self.value_proj = nn.Linear(self.feature_size, self.feature_size)
        self.output_proj = nn.Linear(self.feature_size, self.feature_size)
        self.query_proj.weight.data.normal_(0, 0.01)
        self.key_proj.weight.data.normal_(0, 0.01)
        self.value_proj.weight.data.normal_(0, 0.01)
        self.output_proj.weight.data.normal_(0, 0.01)
        
    def forward(self, input1, input2):
        """
        Args:
            input1: Tensor of shape (batch, channels, feature_size)
            input2: Tensor of shape (batch, channels, feature_size)
        
        Returns:
            Tensor of shape (batch, channels, feature_size), attended features.
        """
        # Project inputs
        query = self.query_proj(input1)  # (batch, channels, feature_size)
        key = self.key_proj(input2)     # (batch, channels, feature_size)
        value = self.value_proj(input2) # (batch, channels, feature_size)
        
        # Compute similarity (scaled dot product)
        scores = torch.einsum('bci,bcj->bcij', query, key)  # (batch, channels, feature_size, feature_size)
        scores = scores / (self.feature_size ** 0.5)
        
        # Compute attention weights
        attention_weights = F.softmax(scores, dim=-1)  # (batch, channels, feature_size, feature_size)
        
        # Compute attended values
        attended = torch.einsum('bcij,bcj->bci', attention_weights, value)  # (batch, channels, feature_size)
        
        # Combine with input1
        output = self.output_proj(attended + input1)  # (batch, channels, feature_size)
        
        return output


class DELayer(nn.Module):
    def __init__(self, n_channels, feature_size, kernel_size=3):
        super().__init__()
        num_groups = 4
        self.conv0 = nn.Conv1d(n_channels, 2 * n_channels, kernel_size=1, padding='same', bias=False)
        self.conv1 = nn.Conv1d(2 * n_channels, 2 * n_channels, kernel_size=kernel_size, padding='same', bias=False)
        self.conv2 = nn.Conv1d(2 * n_channels, 2 * n_channels, kernel_size=kernel_size, padding='same', bias=False)
        self.norm1 = nn.GroupNorm(num_groups, 4 * n_channels)
        self.norm2 = nn.GroupNorm(num_groups, 4 * n_channels)
        # self.norm3 = nn.GroupNorm(num_groups, 2 * n_channels)
        # self.scaled_leaky_relu = ScaledLeakyReLU(negative_slope=0.1, scaled_slope=0.9)
        self.attention = DualInputAttention(4 * n_channels, feature_size)
        self.conv0.weight.data.normal_(0, 0.01)
        self.conv1.weight.data.normal_(0, 0.01)
        self.conv2.weight.data.normal_(0, 0.01)

    def forward(self, z0, x):
        # z = self.conv2(self.norm1(F.relu(self.conv1(z0))))
        # out = self.attention (self.norm2(z),self.norm3(x))
        # z = self.norm4(F.relu(out))
        x = F.relu(self.conv0(x))
        x_conv1 = F.relu(self.conv1(x))
        x_conv2 = F.relu(self.conv2(x_conv1))
        x_concat = torch.cat([x_conv1, x_conv2], dim=1)
        x = self.norm1(x_concat)
        # z = self.norm2(z0)
        z = self.norm2(F.relu(self.attention(z0, x)))
        return z

# Define the Bayesian Neural Network with Dropout 
class BNN(nn.Module):
    def __init__(self, initial_channels, hidden_dim, dropout):
        super(BNN, self).__init__()
        self.dropout_value = dropout
        self.fc1 = nn.Linear(4 * initial_channels, 8* hidden_dim)
        self.dropout1 = nn.Dropout(self.dropout_value)  # Dropout for Monte Carlo Dropout
        self.fc2 = nn.Linear(8 * hidden_dim, 4 * hidden_dim)
        self.dropout2 = nn.Dropout(self.dropout_value)
        self.fc3 = nn.Linear(4 * hidden_dim, 1)

    def forward(self, x, training=True):
        x = torch.relu(self.fc1(x))
        x = self.dropout1(x) if training else x   # Apply dropout during training and inference
        x = torch.relu(self.fc2(x))
        x = self.dropout2(x) if training else x   # Apply dropout during training and inference
        return self.fc3(x)

class RegressorModel(nn.Module):
    def __init__(self, initial_channels, feature_size, hidden_dim, dropout):
        super(RegressorModel, self).__init__()

        # Layers
        # self.conv =nn.Conv1d(initial_channels, 2 * initial_channels, kernel_size=1, bias=False)
        # self.norm = nn.GroupNorm(4, 2 * initial_channels)
        self.f = DELayer(n_channels= initial_channels, feature_size=feature_size )
        self.DEQ = DEQFixedPoint(self.f, anderson, tol=1e-4, max_iter=200, beta=1.0)
        self.pool = nn.AdaptiveAvgPool1d(1)
        self.feedforward = BNN(initial_channels, hidden_dim, dropout)

    def forward(self, x, training=True):
        # x = self.norm(F.relu(self.conv(x)))
        x = self.DEQ(x)
        x = self.pool(x[0]).squeeze(-1)             # Apply Deep Equilibrium Model
        # x = torch.flatten(x, start_dim=1)
        output = self.feedforward(x, training)            # Apply feedforward Model for regression
  
        return output

    
learning_rate = 1e-3 
init_channels = window
inner_channels = window
n_hidden_units =  window
set_dropout = 0.4

model = RegressorModel(initial_channels=init_channels, feature_size=n_features, hidden_dim=n_hidden_units, dropout=set_dropout ).to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate) 
scheduler = lr_scheduler.StepLR(optimizer, step_size=10, gamma=0.5)
In [478]:
X= torch.randn((64, window, n_features), device=device)
f = DELayer(n_channels=init_channels, feature_size= n_features).to(device)
DEQ = DEQFixedPoint(f, anderson, tol=1e-4, max_iter=100, beta=0.9).to(device)
out = DEQ(X)[0]
print(out.shape)
out2 = nn.AdaptiveAvgPool1d(1)(out).squeeze(-1)
print(out2.shape)
out3 = model(X)[0]
print(out3.shape)
torch.Size([64, 60, 18])
torch.Size([64, 60])
torch.Size([1])

📌 Validation and Test Functions

In [479]:
# def validation(loss_fn, num_samples=100):
#     model.train()  # Enable stochastic behavior (dropout during inference)
#     total_loss = 0.0
#     total_samples = 0

#     with torch.no_grad():
#         for X, y in valloader:
#             X, y = X.to(device).float(), y.to(device).float()

#             stochastic_predictions = []

#             # Perform stochastic forward passes
#             for _ in range(num_samples):
#                 y_pred = model(X, training=True).squeeze()
#                 stochastic_predictions.append(y_pred.cpu().numpy())

#             # Convert to NumPy array
#             stochastic_predictions = np.array(stochastic_predictions)  # Shape: (num_samples, batch_size)

#             # Calculate mean prediction
#             mean_pred = torch.tensor(np.mean(stochastic_predictions, axis=0)).to(device)

#             # Compute batch loss
#             loss = loss_fn(mean_pred, y)
#             total_loss += loss.item() * X.size(0)  # Accumulate the weighted batch loss
#             total_samples += X.size(0)

#     avg_val_loss = total_loss / total_samples  # Calculate average loss
#     return avg_val_loss

def validation(loss_fn):
    model.eval()  # Set the model to evaluation mode
    total_loss = 0.0
    total_samples = 0

    with torch.no_grad():  # Disable gradient computation for validation
        for X, y in valloader:  # Iterate over all batches in the validation loader
            X, y = X.to(device).to(torch.float32), y.to(device).to(torch.float32)
            
            # Use enable_grad for specific models if required
            if isinstance(model, RegressorModel):  # Replace with the DEQ model class
                with torch.enable_grad():
                    y_pred = model(X).squeeze()
            else:
                y_pred = model(X).squeeze()
            
            y = y.squeeze()
            loss = loss_fn(y_pred, y)
            total_loss += loss.item() * X.size(0)  # Accumulate the weighted batch loss
            total_samples += X.size(0)  # Count the total samples processed

    avg_val_loss = total_loss / total_samples  # Calculate average loss
    return avg_val_loss

loss_L1 = nn.L1Loss()
    
def test():
    model.train()  # Enable stochastic behavior (dropout during inference)
    
    total_loss_MSE = 0.0
    total_loss_L1 = 0.0
    total_ASUE = 0.0
    all_stochastic_predictions = []
    all_pred_values = []
    all_true_values = []
    all_uncertainties = []  # Collect uncertainties
    num_samples = 100

    with torch.no_grad():
        for X, y in testloader:
            X, y = X.to(device).float(), y.to(device).float()

            stochastic_predictions = []

            # Perform stochastic forward passes
            for _ in range(num_samples):
                y_pred = model(X, training=True).squeeze()
                stochastic_predictions.append(y_pred.cpu().numpy())
            
            # Convert to NumPy array
            stochastic_predictions = np.array(stochastic_predictions)  # Shape: (num_samples, batch_size)

            # Calculate mean and uncertainty
            mean_pred = torch.tensor(np.mean(stochastic_predictions, axis=0)).to(device)
            uncertainty_batch = np.var(stochastic_predictions, axis=0)  # Variance as uncertainty

            # Compute batch losses
            loss_MSE = torch.mean((mean_pred - y) ** 2).item()
            loss_L1_val = loss_L1(mean_pred, y).item()
            ASUE = torch.mean(torch.relu(y - mean_pred)).item()

            # Aggregate results
            total_loss_MSE += loss_MSE * X.size(0)
            total_loss_L1 += loss_L1_val * X.size(0)
            total_ASUE += ASUE * X.size(0)

            # Collect predictions, true values, and uncertainties
            all_stochastic_predictions.append(stochastic_predictions)
            all_pred_values.append(mean_pred.cpu().numpy())
            all_true_values.append(y.cpu().numpy())
            all_uncertainties.append(uncertainty_batch)
        
    # Average metrics
    total_samples = len(testloader.dataset)
    avg_loss_MSE = total_loss_MSE / total_samples
    avg_loss_L1 = total_loss_L1 / total_samples
    avg_ASUE = total_ASUE / total_samples

    # Concatenate all predictions, uncertainties, and true values
    all_pred_values = np.concatenate([x.flatten() for x in all_pred_values], axis=0)
    all_stochastic_predictions = np.concatenate(all_stochastic_predictions, axis=1)
    all_true_values = np.concatenate(all_true_values, axis=0)
    all_uncertainties = np.concatenate(all_uncertainties, axis=0)  # Ensure consistency

    return avg_loss_MSE, avg_loss_L1, avg_ASUE, all_uncertainties, all_stochastic_predictions, all_pred_values, all_true_values

📌 Training Loop: I have trained using Adam Optimizer for 35 epochs with learning rate = 0.01 and scheduler with degration 0.5 after 8 epochs

In [ ]:
# Initialize lists for tracking losses
T, V = [], []
epochs = 35

loss_fn = nn.MSELoss()
siterations, sresiduals = [], []
stats = {}

# Define forward hook for DEQ statistics
def forward_hook(module, input, output):
    stats[module] = output

# Register forward hook
hook_handle = model.DEQ.register_forward_hook(forward_hook)

best_val_loss = float('inf')  # Track best validation loss

pbar = tqdm(range(epochs), desc="Training Progress", dynamic_ncols=True)
for epoch in pbar:
    model.train()
    epoch_loss = 0
    sresiduals_list, siterations_list = [], []

    for X, y in trainloader:
        X, y = X.to(device, dtype=torch.float32), y.to(device, dtype=torch.float32)

        # Forward pass
        y_pred = model(X).squeeze()
        y = y.squeeze()
        loss = loss_fn(y_pred, y)
        epoch_loss += loss.item()

        # Backward pass
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

        # Collect DEQ statistics  
        iterations = stats.get(model.DEQ, (None, [0, 0]))[1] 
        if iterations and len(iterations) >= 2:
            sresiduals_list.append(iterations[0][-1])
            siterations_list.append(iterations[1])

    # Step the scheduler
    scheduler.step()

    # Validation loss
    val_loss = validation(loss_fn)
    T.append(epoch_loss / len(trainloader))
    V.append(val_loss)
 
    # Update tqdm bar
    pbar.set_postfix({'Train Loss': T[-1], 'Val Loss': val_loss})

    # Track DEQ statistics
    sresiduals.append(np.mean(sresiduals_list) if sresiduals_list else 0)
    siterations.append(np.mean(siterations_list) if siterations_list else 0)

    # Save best model
    if val_loss < best_val_loss:
        best_val_loss = val_loss
        torch.save(model.state_dict(), 'best_model.pth')

    print(f'Epoch {epoch + 1}/{epochs} - Train Loss: {T[-1]:.4f}, Val Loss: {V[-1]:.4f}')

# Clean up hook handle
hook_handle.remove()
In [ ]:
print(sresiduals)
print(siterations)
In [ ]:
torch.save(model.state_dict(), f'saved_models\\DEM_{file_name}')
In [ ]:
# Perform testing with Monte Carlo Dropout
mse, l1, asue, uncertainty, all_y_pred, y_pred, y = test()

# Calculate RMSE and scale by RUL max value
rmse = rul_max * np.sqrt(mse)

# Print results
print(f"Test Results:\n"
      f"RMSE: {round(rmse, 2)}\n"
      f"L1 Loss: {round(l1, 2)}\n"
      f"ASUE: {round(asue, 2)}")

# Mean uncertainty
mean_uncertainty = np.mean(uncertainty)
print(f"Mean Uncertainty: {round(mean_uncertainty, 4)}")

# Optional: Visualize predictions, true values, and uncertainties
import matplotlib.pyplot as plt

plt.figure(figsize=(14, 8))

# True vs predicted RUL
plt.plot(range(len(y)), y, label="True Values", c="lightseagreen", marker=".", alpha=0.7)
plt.plot(range(len(y_pred)), y_pred, label="Predicted Values", c="salmon", marker=".", alpha=0.7)

# Uncertainty as error bars
plt.fill_between(
    range(len(y_pred)),
    y_pred - np.sqrt(uncertainty),
    y_pred + np.sqrt(uncertainty),
    color="salmon", alpha=0.2, label="Uncertainty"
)

# Add labels and legend
plt.xlabel("Sample Index", fontsize=14)
plt.ylabel("RUL", fontsize=14)
plt.title("True vs Predicted RUL with Uncertainty", fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.savefig(f'figures\\Responses_with_Uncertainty_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
In [ ]:
print(y)
print(y_pred)
print(1*(y_pred>y))
print(all_y_pred.shape)
In [ ]:
print(all_y_pred)
In [ ]:
# Create a figure and axes
fig, ax = plt.subplots(1, 1, figsize=(14, 8))

# Plot predictions
ax.plot(np.arange(1, y_pred.shape[0] + 1), y_pred, label='Predictions', color='salmon', marker='.')

# Plot true values
ax.plot(np.arange(1, y.shape[0] + 1), y, label='True Values', color='lightseagreen', marker='.')

# Set limits, labels, grid, and legend
ax.set_ylim([0, 1])
ax.set_xlabel('Test Engine Units', fontsize=16)
ax.set_ylabel('RUL', fontsize=16)
ax.grid(True)
ax.legend()
plt.savefig(f'figures\\Responses_without_Uncertainty_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
# Show the plot
plt.show()
In [ ]:
plt.plot(np.arange(1,len(T)+1), T, label= 'Train loss')
plt.plot(np.arange(1,len(V)+1), V, label = 'Validation loss')
plt.legend()
plt.grid()
plt.xlabel('Epochs')
plt.ylabel('MSE')
plt.savefig(f'figures\\Training_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show
In [ ]:
fig, ax = plt.subplots(figsize = (12,8))

def animate(i):
    ax.clear()
    line1, = ax.plot(np.arange(1,i+1), T[:i], label = 'train_loss')
    line2, = ax.plot(np.arange(1,i+1), V[:i], label = 'val_loss')
    ax.legend()
    ax.grid(True)
    ax.set_xlim(0,101)
    ax.set_ylim(0,4000)
    ax.set_xlabel('epochs')
    ax.set_ylabel('MSE')
                     
    return line1, line2
                     
    
animation = FA(fig, animate, np.arange(1,len(T)+1), interval = 50)

%time animation.save('animation3.gif', writer='imagemagick', fps=20)

plt.close(fig)
In [ ]:
# Count total parameters
total_params = sum(p.numel() for p in model.parameters())
print(f"Total number of parameters: {total_params}")

# Count trainable parameters only
trainable_params = sum(p.numel() for p in model.parameters() if p.requires_grad)
print(f"Trainable parameters: {trainable_params}")
In [ ]:
from torchinfo import summary
summary(model, input_size=(64, window, n_features))

3. Prediction on Test Set¶

In [ ]:
model.load_state_dict(torch.load(f'saved_models\\DEM_{file_name}', map_location=torch.device(device)))

Phm08 Score Metric

In [ ]:
def calculate_phm08_score(true_rul, predicted_rul, alpha=13, beta=10):
    """
    Calculate the PHM08 score metric for RUL prediction.

    Parameters:
        true_rul (array-like): Array of true RUL values.
        predicted_rul (array-like): Array of predicted RUL values.
        alpha (float): Scaling factor for early predictions (default=13).
        beta (float): Scaling factor for late predictions (default=10).

    Returns:
        float: Total PHM08 score.
    """
    true_rul = np.array(true_rul)
    predicted_rul = np.array(predicted_rul)

    if true_rul.shape != predicted_rul.shape:
        raise ValueError("Shape mismatch between true and predicted RUL arrays")

    d = (predicted_rul - true_rul) * MAX_RUL

    score = np.where(d < 0, np.exp(-d / alpha) - 1, np.exp(d / beta) - 1)
    return np.sum(score)
In [ ]:
num_samples = 100

def test(num_samples=num_samples):
    model.train()  # Enable stochastic behavior (dropout during inference)
    
    total_loss_MSE = 0.0
    total_loss_L1 = 0.0
    total_ASUE = 0.0
    total_phm08_score = 0.0
    all_stochastic_predictions = []
    all_pred_values = []
    all_true_values = []
    all_uncertainties = []  # Collect uncertainties

    with torch.no_grad():
        for X, y in testloader:
            X, y = X.to(device).float(), y.to(device).float()

            stochastic_predictions = []

            # Perform stochastic forward passes
            for _ in range(num_samples):
                y_pred = model(X, training=True).squeeze()
                stochastic_predictions.append(y_pred.cpu().numpy())
            
            # Convert to NumPy array
            stochastic_predictions = np.array(stochastic_predictions)  # Shape: (num_samples, batch_size)

            # Calculate mean and uncertainty with std
            mean_pred = torch.tensor(np.mean(stochastic_predictions, axis=0)).to(device)
            uncertainty_batch = np.std(stochastic_predictions, axis=0)  # Variance as uncertainty

            # Compute batch losses
            loss_MSE = torch.mean((mean_pred - y) ** 2).item()
            loss_L1_val = loss_L1(mean_pred, y).item()
            ASUE = torch.mean(torch.relu(y - mean_pred)).item()

            # Compute PHM08 score
            phm08_score = calculate_phm08_score(y.cpu().numpy(), mean_pred.cpu().numpy())

            # Aggregate results
            total_loss_MSE += loss_MSE * X.size(0)
            total_loss_L1 += loss_L1_val * X.size(0)
            total_ASUE += ASUE * X.size(0)
            total_phm08_score += phm08_score

            # Collect predictions, true values, and uncertainties
            all_stochastic_predictions.append(stochastic_predictions)
            all_pred_values.append(mean_pred.cpu().numpy())
            all_true_values.append(y.cpu().numpy())
            all_uncertainties.append(uncertainty_batch)
        
    # Average metrics
    total_samples = len(testloader.dataset)
    avg_loss_MSE = total_loss_MSE / total_samples
    avg_loss_L1 = total_loss_L1 / total_samples
    avg_ASUE = total_ASUE / total_samples
    avg_phm08_score = total_phm08_score

    # Concatenate all predictions, uncertainties, and true values
    all_pred_values = np.concatenate([x.flatten() for x in all_pred_values], axis=0)
    all_stochastic_predictions = np.concatenate(all_stochastic_predictions, axis=1)
    all_true_values = np.concatenate(all_true_values, axis=0)
    all_uncertainties = np.concatenate(all_uncertainties, axis=0)  # Ensure consistency

    return avg_loss_MSE, avg_loss_L1, avg_ASUE, avg_phm08_score, all_uncertainties, all_stochastic_predictions, all_pred_values, all_true_values
In [ ]:
# Perform testing with Monte Carlo Dropout
mse, l1, asue, score, uncertainty, all_y_pred, y_pred, y = test()

# Calculate RMSE and scale by RUL max value
rmse = rul_max * np.sqrt(mse)

# Print results
print(f"Test Results:\n"
      f"RMSE: {round(rmse, 2)}\n"
      f"Score: {round(score, 2)}\n"
      f"L1 Loss: {round(l1, 2)}\n"
      f"ASUE: {round(asue, 2)}")

# Mean uncertainty
mean_uncertainty = np.mean(uncertainty)
print(f"Mean Uncertainty: {round(mean_uncertainty, 4)}")

# Optional: Visualize predictions, true values, and uncertainties
import matplotlib.pyplot as plt

plt.figure(figsize=(14, 8))

# True vs predicted RUL
plt.plot(range(len(y)), y, label="True Values", c="lightseagreen", marker=".", alpha=0.7)
plt.plot(range(len(y_pred)), y_pred, label="Predicted Values", c="salmon", marker=".", alpha=0.7)

# Uncertainty as error bars
plt.fill_between(
    range(len(y_pred)),
    y_pred - uncertainty,
    y_pred + uncertainty,
    color="salmon", alpha=0.2, label="Uncertainty"
)

# Add labels and legend
plt.xlabel("Sample Index", fontsize=14)
plt.ylabel("RUL", fontsize=14)
plt.title("True vs Predicted RUL with Uncertainty", fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.savefig(f'figures\\Responses_with_Uncertainty_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
In [ ]:
values = np.array([y, y_pred, uncertainty]).T
df_values = pd.DataFrame(values)
values = df_values.sort_values(0, ascending=False).to_numpy()

y, y_pred, uncertainty = np.array(values).T

# y, y_pred, uncertainty = MAX_RUL * np.array(values).T
# uncertainty = MAX_RUL * uncertainty
In [ ]:
plt.hist(uncertainty, bins=50, color="blue", alpha=0.7, label="Uncertainty Distribution")
plt.xlabel("Uncertainty")
plt.ylabel("Frequency")
plt.title("Histogram of Prediction Uncertainty")
plt.legend()
plt.grid(True)
plt.show()
In [ ]:
# confidence = 0.99  # Confidence level
# zeta = 2.576

# # Calculate confidence interval
# standard_error = uncertainty / np.sqrt(num_samples)  # Standard error of the mean
# margin_of_error = zeta * standard_error  # Margin of error

confidence = 0.95  # Confidence level
margin_of_error = confidence * uncertainty  # Margin of error
In [ ]:
plt.figure(figsize=(8, 6)) # Adjust figure size for publication

# True Values (more prominent line)
plt.plot(range(len(y)), y, label="True RUL", c="lightseagreen", linewidth=2)

# Predicted Values (smaller scatter points)
plt.scatter(range(len(y_pred)), y_pred, label="Predicted RUL", c="salmon", edgecolor="black", alpha=0.7, zorder=3, s=15) # reduced scatter point size to 15

# Vertical Lines (improved clarity)
for i in range(len(y)):
    plt.plot([i, i], [y_pred[i], y[i]], c="gray", linestyle="--", linewidth=0.8, alpha=0.7)

# Confidence Interval (adjusted opacity, more descriptive label)
plt.fill_between(
    range(len(y_pred)),
    y_pred - margin_of_error,
    y_pred + margin_of_error,
    color="salmon",
    alpha=0.3,  # Reduced opacity
    label=f"{int(confidence * 100)}% Prediction Interval"  # More accurate label
)

# Labels and Title (improved clarity)
dataset_name = file_name[:-4]
plt.xlabel("Testing Sample Index") # Removed unnecessary context from label
plt.ylabel("Remaining Useful Life (RUL)") # Removed unnecessary context from label
plt.title(f"Comparison of True and Predicted RUL for {dataset_name}") # Removed unnecessary context from title

# Legend and Grid
plt.legend()

# Save and Show
plt.savefig(f'figures/True_vs_PredictedRUL_{file_name}.png')
plt.show()
In [ ]:
# Calculate errors and mean predictions
errors = np.abs(y_pred - y)

# Scatter plot with enhancements
plt.figure(figsize=(8, 6))
sc = plt.scatter(uncertainty, errors, c=errors, cmap='viridis', alpha=0.7, edgecolor='k')

# Add colorbar
cbar = plt.colorbar(sc)
cbar.set_label('Error Magnitude')

# Add grid and labels
plt.grid(alpha=0.3)
plt.xlabel('Uncertainty', fontsize=12)
plt.ylabel('Error', fontsize=12)
plt.title('Uncertainty vs. Error', fontsize=14)

# Polynomial fit (degree 2 or higher)
degree = 1  # try 3, 4, etc. for more complex trends
z = np.polyfit(uncertainty, errors, degree)
p = np.poly1d(z)

# Sort for smooth plotting
sorted_idx = np.argsort(uncertainty)
plt.plot(uncertainty[sorted_idx], p(uncertainty[sorted_idx]),
         color='red', linestyle='--', label=f'Polynomial Regression (deg {degree})')

# plt.plot(uncertainty, p(uncertainty), color='red', linestyle='--', label='Trend Line')

# Add legend
plt.legend()

# Show the plot
plt.tight_layout()
plt.savefig(f'figures\\UncertaintyVSError_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()


# Compute statistics
correlation = np.corrcoef(uncertainty, errors)[0, 1]
print(f"Correlation between uncertainty and error: {correlation:.2f}")

# Evaluate performance in high/low uncertainty regions
threshold = np.percentile(uncertainty, 75)  # Upper 25% uncertainty
high_uncertainty_mask = uncertainty > threshold
low_uncertainty_mask = uncertainty <= threshold

high_uncertainty_error = errors[high_uncertainty_mask].mean()
low_uncertainty_error = errors[low_uncertainty_mask].mean()

print(f"Mean error in high uncertainty regions: {high_uncertainty_error:.2f}")
print(f"Mean error in low uncertainty regions: {low_uncertainty_error:.2f}")
In [ ]:
# Calculate absolute error
errors = np.abs(y_pred - y)

# Plot
plt.figure(figsize=(8, 6))
sc = plt.scatter(y, errors, c=errors, cmap='magma', alpha=0.7, edgecolor='k')

# Add colorbar
cbar = plt.colorbar(sc)
cbar.set_label('Error Magnitude')

# Add labels and title
plt.xlabel('True Values', fontsize=12)
plt.ylabel('Absolute Error', fontsize=12)
plt.title('True Values vs. Prediction Error', fontsize=14)
plt.grid(alpha=0.3)

# Optional trend line
z = np.polyfit(y, errors, 2)
p = np.poly1d(z)
sorted_idx = np.argsort(y)
plt.plot(y[sorted_idx], p(y[sorted_idx]), color='red', linestyle='--', label='Trend (deg 2)')

# Add legend
plt.legend()

# Save and show
plt.tight_layout()
plt.savefig(f'figures\\TrueValueVsError_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()

4. Prediction on Test Set Without Dropout in the Inference¶

In [ ]:
model.load_state_dict(torch.load(f'saved_models\\DEM_{file_name}', map_location=torch.device(device)))
In [ ]:
def calculate_phm08_score(true_rul, predicted_rul, alpha=13, beta=10):
    """
    Calculate the PHM08 score metric for RUL prediction.

    Parameters:
        true_rul (array-like): Array of true RUL values.
        predicted_rul (array-like): Array of predicted RUL values.
        alpha (float): Scaling factor for early predictions (default=13).
        beta (float): Scaling factor for late predictions (default=10).

    Returns:
        float: Total PHM08 score.
    """
    true_rul = np.array(true_rul)
    predicted_rul = np.array(predicted_rul)

    if true_rul.shape != predicted_rul.shape:
        raise ValueError("Shape mismatch between true and predicted RUL arrays")

    d = (predicted_rul - true_rul) * MAX_RUL

    score = np.where(d < 0, np.exp(-d / alpha) - 1, np.exp(d / beta) - 1)
    return np.sum(score)
In [ ]:
def test_without_dropout():
    model.eval()  # Disable stochastic behavior (dropout during inference)

    total_loss_MSE = 0.0
    total_loss_L1 = 0.0
    total_ASUE = 0.0
    total_phm08_score = 0.0
    all_pred_values = []
    all_true_values = []

    with torch.no_grad():
        for X, y in testloader:
            X, y = X.to(device).float(), y.to(device).float()

            # Perform a single deterministic forward pass
            y_pred = model(X).squeeze()

            # Compute batch losses
            loss_MSE = torch.mean((y_pred - y) ** 2).item()
            loss_L1_val = loss_L1(y_pred, y).item()
            ASUE = torch.mean(torch.relu(y - y_pred)).item()

            # Compute PHM08 score
            phm08_score = calculate_phm08_score(y.cpu().numpy(), y_pred.cpu().numpy())

            # Aggregate results
            total_loss_MSE += loss_MSE * X.size(0)
            total_loss_L1 += loss_L1_val * X.size(0)
            total_ASUE += ASUE * X.size(0)
            total_phm08_score += phm08_score

            # Collect predictions and true values
            all_pred_values.append(y_pred.cpu().numpy())
            all_true_values.append(y.cpu().numpy())

    # Average metrics
    total_samples = len(testloader.dataset)
    avg_loss_MSE = total_loss_MSE / total_samples
    avg_loss_L1 = total_loss_L1 / total_samples
    avg_ASUE = total_ASUE / total_samples
    avg_phm08_score = total_phm08_score

    # Concatenate all predictions and true values
    all_pred_values = np.concatenate([x.flatten() for x in all_pred_values], axis=0)
    all_true_values = np.concatenate(all_true_values, axis=0)

    return avg_loss_MSE, avg_loss_L1, avg_ASUE, avg_phm08_score, all_pred_values, all_true_values
In [ ]:
# Perform testing without Monte Carlo Dropout
mse, l1, asue, score, y_pred, y = test_without_dropout()

# Calculate RMSE and scale by RUL max value
rmse = rul_max * np.sqrt(mse)

# Print results
print(f"Test Results:\n"
      f"RMSE: {rmse:.2f}\n"
      f"Score: {score:.2f}\n"
      f"L1 Loss: {l1:.2f}\n"
      f"ASUE: {asue:.3f}")

5. Methods for Weighting Predictions by Uncertainty¶

In [ ]:
# Perform testing with Monte Carlo Dropout
mse, l1, asue, score, uncertainty, all_y_pred, y_pred, y = test()
In [ ]:
loss_MSE = np.mean((y_pred - y) ** 2).item()
rmse = rul_max * np.sqrt(mse)
print(f'RMSE:{rmse:.2f}')
In [ ]:
all_y_pred_mean = all_y_pred.mean(axis=0)
weights = abs(all_y_pred - all_y_pred_mean)
weights = np.exp(-weights)
# weights = 1 / weights
weights /= weights.sum(axis=0)
weighted_predictions = np.sum(weights * all_y_pred, axis=0)
In [ ]:
mse = np.mean((weighted_predictions - y) ** 2).item()
rmse = rul_max * np.sqrt(mse)
print(f'RMSE:{rmse:.2f}')

Threshold In Uncertainty

In [ ]:
# Define a threshold based on deviations (e.g., 1.5x the mean deviation)
threshold = 1.5 * weights.mean(axis=0)

# Mask predictions with deviations exceeding the threshold
valid_predictions = np.where(weights <= threshold, all_y_pred, np.nan)

# Compute predictions as the mean of valid predictions
weighted_predictions = np.nanmean(valid_predictions, axis=0)

mse = np.mean((weighted_predictions - y) ** 2).item()
rmse = rul_max * np.sqrt(mse)
print(f'RMSE:{rmse:.2f}')

4. Temperature Scaling - Calibration on Uncertaintity¶

In [480]:
model.load_state_dict(torch.load(f'saved_models\\DEM_{file_name}', map_location=torch.device(device)))
C:\Users\spiro\AppData\Local\Temp\ipykernel_22120\4004060783.py:1: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature.
  model.load_state_dict(torch.load(f'saved_models\\DEM_{file_name}', map_location=torch.device(device)))
Out[480]:
<All keys matched successfully>
In [481]:
def calculate_phm08_score(true_rul, predicted_rul, alpha=13, beta=10):
    """
    Calculate the PHM08 score metric for RUL prediction.

    Parameters:
        true_rul (array-like): Array of true RUL values.
        predicted_rul (array-like): Array of predicted RUL values.
        alpha (float): Scaling factor for early predictions (default=13).
        beta (float): Scaling factor for late predictions (default=10).

    Returns:
        float: Total PHM08 score.
    """
    true_rul = np.array(true_rul)
    predicted_rul = np.array(predicted_rul)

    if true_rul.shape != predicted_rul.shape:
        raise ValueError("Shape mismatch between true and predicted RUL arrays")

    d = (predicted_rul - true_rul) * MAX_RUL

    score = np.where(d < 0, np.exp(-d / alpha) - 1, np.exp(d / beta) - 1)
    return np.sum(score)
In [482]:
num_samples = 100

def test(num_samples=num_samples):
    model.train()  # Enable stochastic behavior (dropout during inference)
    
    total_loss_MSE = 0.0
    total_loss_L1 = 0.0
    total_ASUE = 0.0
    total_phm08_score = 0.0
    all_stochastic_predictions = []
    all_pred_values = []
    all_true_values = []
    all_uncertainties = []  # Collect uncertainties

    with torch.no_grad():
        for X, y in testloader:
            X, y = X.to(device).float(), y.to(device).float()

            stochastic_predictions = []

            # Perform stochastic forward passes
            for _ in range(num_samples):
                y_pred = model(X, training=True).squeeze()
                stochastic_predictions.append(y_pred.cpu().numpy())
            
            # Convert to NumPy array
            stochastic_predictions = np.array(stochastic_predictions)  # Shape: (num_samples, batch_size)

            # Calculate mean and uncertainty with std
            mean_pred = torch.tensor(np.mean(stochastic_predictions, axis=0)).to(device)
            uncertainty_batch = np.std(stochastic_predictions, axis=0)  # Variance as uncertainty

            # Compute batch losses
            loss_MSE = torch.mean((mean_pred - y) ** 2).item()
            loss_L1_val = loss_L1(mean_pred, y).item()
            ASUE = torch.mean(torch.relu(y - mean_pred)).item()

            # Compute PHM08 score
            phm08_score = calculate_phm08_score(y.cpu().numpy(), mean_pred.cpu().numpy())

            # Aggregate results
            total_loss_MSE += loss_MSE * X.size(0)
            total_loss_L1 += loss_L1_val * X.size(0)
            total_ASUE += ASUE * X.size(0)
            total_phm08_score += phm08_score

            # Collect predictions, true values, and uncertainties
            all_stochastic_predictions.append(stochastic_predictions)
            all_pred_values.append(mean_pred.cpu().numpy())
            all_true_values.append(y.cpu().numpy())
            all_uncertainties.append(uncertainty_batch)
        
    # Average metrics
    total_samples = len(testloader.dataset)
    avg_loss_MSE = total_loss_MSE / total_samples
    avg_loss_L1 = total_loss_L1 / total_samples
    avg_ASUE = total_ASUE / total_samples
    avg_phm08_score = total_phm08_score

    # Concatenate all predictions, uncertainties, and true values
    all_pred_values = np.concatenate([x.flatten() for x in all_pred_values], axis=0)
    all_stochastic_predictions = np.concatenate(all_stochastic_predictions, axis=1)
    all_true_values = np.concatenate(all_true_values, axis=0)
    all_uncertainties = np.concatenate(all_uncertainties, axis=0)  # Ensure consistency

    return avg_loss_MSE, avg_loss_L1, avg_ASUE, avg_phm08_score, all_uncertainties, all_stochastic_predictions, all_pred_values, all_true_values
In [483]:
def validate_4_calibrate(num_samples=100):
    model.train()  # Dropout needs to stay on
    
    all_pred_values = []
    all_uncertainties = []
    all_true_values = []

    with torch.no_grad():
        for X, y in valloader:
            X, y = X.to(device).float(), y.to(device).float()

            stochastic_predictions = []

            for _ in range(num_samples):
                y_pred = model(X, training=True).squeeze()
                stochastic_predictions.append(y_pred.cpu().numpy())

            stochastic_predictions = np.array(stochastic_predictions)  # (num_samples, batch_size)
            mean_pred = np.mean(stochastic_predictions, axis=0)
            uncertainty = np.std(stochastic_predictions, axis=0)

            all_pred_values.append(mean_pred)
            all_uncertainties.append(uncertainty)
            all_true_values.append(y.cpu().numpy())

    # Flatten
    val_pred_means = np.concatenate(all_pred_values)
    val_pred_stds = np.concatenate(all_uncertainties)
    val_targets = np.concatenate(all_true_values)

    return val_pred_means, val_pred_stds, val_targets



val_pred_means, val_pred_stds, val_targets = validate_4_calibrate(num_samples=100)
In [484]:
mse, l1, asue, score, uncertainty, all_y_pred, y_pred, y = test()
In [485]:
def compute_coverage(y_true, y_pred, scaled_stds):
    """
    Computes the empirical coverage: the fraction of samples where
    the true value is inside the predicted confidence interval.
    
    Parameters:
    - y_true : np.array, ground truth values
    - y_pred : np.array, predicted means
    - scaled_stds : np.array, calibrated standard deviations (after applying k-scaling)

    Returns:
    - coverage : float, percentage of samples inside the confidence interval
    """
    lower_bounds = y_pred - scaled_stds
    upper_bounds = y_pred + scaled_stds

    # Check if the true value lies inside the interval
    inside = (y_true >= lower_bounds) & (y_true <= upper_bounds)

    # Calculate coverage as a percentage
    coverage = np.mean(inside) * 100.0
    return coverage


coverage = compute_coverage(y, y_pred, uncertainty)
print(f"Ratio of coverage:{coverage}%")
Ratio of coverage:63.0%
In [486]:
values = np.array([y, y_pred, uncertainty]).T
df_values = pd.DataFrame(values)
values = df_values.sort_values(0, ascending=False).to_numpy()

y, y_pred, uncertainty = np.array(values).T

# y, y_pred, uncertainty = MAX_RUL * np.array(values).T
# uncertainty = MAX_RUL * uncertainty

Simple scaling-based calibration method

In [487]:
import scipy.stats as st

import scipy.optimize as optim

def calibrate_scaling(val_pred_means, val_stds, val_targets, alpha=0.05):
    def coverage_diff(k):
        lower = val_pred_means - k * val_stds
        upper = val_pred_means + k * val_stds
        coverage = np.mean((val_targets >= lower) & (val_targets <= upper))
        return abs((1 - alpha) - coverage)

    result = optim.minimize_scalar(coverage_diff, bounds=(0.1, 5.0), method='bounded')  # Adjust bounds
    return result.x

k_opt = calibrate_scaling(val_pred_means, val_pred_stds, val_targets, alpha=0.05)
display(k_opt)
2.6742587996512373
In [488]:
calibrated_lower = y_pred - k_opt * uncertainty
calibrated_upper = y_pred + k_opt * uncertainty

test_coverage = np.mean((y >= calibrated_lower) & (y <= calibrated_upper))
display(f'Coverage {test_coverage*100}%')
'Coverage 94.0%'
In [489]:
def plot_reliability(pred_means, pred_stds, targets):
    confidences = np.linspace(0.5, 0.99, 20)
    coverages = []

    for conf in confidences:
        z = st.norm.ppf(1 - (1 - conf)/2)
        lower = pred_means - z * pred_stds
        upper = pred_means + z * pred_stds
        coverage = np.mean((targets >= lower) & (targets <= upper))
        coverages.append(coverage)

    plt.plot(confidences, coverages, label='Empirical Coverage')
    plt.plot([0.5, 0.99], [0.5, 0.99], 'k--', label='Ideal Calibration')
    plt.xlabel("Nominal Confidence")
    plt.ylabel("Empirical Coverage")
    plt.legend()
    plt.grid()
    plt.title("Calibration Curve")
    plt.show()

plot_reliability(y_pred, uncertainty, y)
No description has been provided for this image
In [490]:
# confidence = 0.99  # Confidence level
# zeta = 2.576

# # Calculate confidence interval
# standard_error = uncertainty / np.sqrt(num_samples)  # Standard error of the mean
# margin_of_error = zeta * standard_error  # Margin of error

confidence = 0.95  # Confidence level
margin_of_error = k_opt * uncertainty  # Margin of error
margin_of_error = confidence * margin_of_error   # Margin of error
In [491]:
# Boolean mask: True if point is inside the confidence interval
inside_interval = (y >= y_pred - margin_of_error) & (y <= y_pred + margin_of_error)

# Assign colors: green if inside, red if outside
colors_interval = np.where(inside_interval, 'mediumseagreen', 'crimson')

plt.figure(figsize=(8, 6)) # Adjust figure size for publication
plt.ylim(0, MAX_RUL+1)

# True Values (more prominent line)
plt.plot(range(len(y)), y * MAX_RUL, label="True RUL", c="lightseagreen", linewidth=2)

# Predicted Values (smaller scatter points)
plt.scatter(range(len(y_pred)), y_pred * MAX_RUL, label="Predicted RUL", c=colors_interval, edgecolor="black", alpha=0.7, zorder=3, s=15) # reduced scatter point size to 15

# Vertical Lines (improved clarity)
for i in range(len(y)):
    plt.plot([i, i], [y_pred[i], y[i]], c="gray", linestyle="--", linewidth=0.8, alpha=0.7)

# Confidence Interval (adjusted opacity, more descriptive label)
plt.fill_between(
    range(len(y_pred)),
    (y_pred - margin_of_error)* MAX_RUL,
    np.minimum(y_pred + margin_of_error, 1) * MAX_RUL,
    color="salmon",
    alpha=0.3,  # Reduced opacity
    label=f"{int(confidence * 100)}% Prediction Interval"  # More accurate label
)

# Labels and Title (improved clarity)
dataset_name = file_name[:-4]
plt.xlabel("Testing Sample Index") # Removed unnecessary context from label
plt.ylabel("Remaining Useful Life (RUL)") # Removed unnecessary context from label
plt.title(f"Comparison of True and Predicted RUL for {dataset_name}") # Removed unnecessary context from title

# Legend and Grid
plt.legend()

# Save and Show
plt.savefig(f'figures/True_vs_PredictedRUL_simple_calibration_{file_name}.png')
plt.show()
No description has been provided for this image
In [492]:
# Calculate errors and mean predictions
errors = np.abs(y_pred - y)

# Scatter plot with enhancements
plt.figure(figsize=(8, 6))
sc = plt.scatter(margin_of_error, errors, c=errors, cmap='viridis', alpha=0.7, edgecolor='k')

# Add colorbar
cbar = plt.colorbar(sc)
cbar.set_label('Error Magnitude')

# Add grid and labels
plt.grid(alpha=0.3)
plt.xlabel('Uncertainty', fontsize=12)
plt.ylabel('Error', fontsize=12)
plt.title('Uncertainty vs. Error', fontsize=14)

# Polynomial fit (degree 2 or higher)
degree = 1  # try 3, 4, etc. for more complex trends
z = np.polyfit(margin_of_error, errors, degree)
p = np.poly1d(z)

# Sort for smooth plotting
sorted_idx = np.argsort(margin_of_error)
plt.plot(margin_of_error[sorted_idx], p(margin_of_error[sorted_idx]),
         color='red', linestyle='--', label=f'Polynomial Regression (deg {degree})')

# plt.plot(uncertainty, p(uncertainty), color='red', linestyle='--', label='Trend Line')

# Add legend
plt.legend()

# Show the plot
plt.tight_layout()
plt.savefig(f'figures\\UncertaintyVSError_simple_calibration_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()


# Compute statistics
correlation = np.corrcoef(margin_of_error, errors)[0, 1]
print(f"Correlation between uncertainty and error: {correlation:.2f}")

# Evaluate performance in high/low uncertainty regions
threshold = np.percentile(margin_of_error, 75)  # Upper 25% uncertainty
high_uncertainty_mask = margin_of_error > threshold
low_uncertainty_mask = margin_of_error <= threshold

high_uncertainty_error = errors[high_uncertainty_mask].mean()
low_uncertainty_error = errors[low_uncertainty_mask].mean()

print(f"Mean error in high uncertainty regions: {high_uncertainty_error:.2f}")
print(f"Mean error in low uncertainty regions: {low_uncertainty_error:.2f}")
No description has been provided for this image
Correlation between uncertainty and error: 0.22
Mean error in high uncertainty regions: 0.06
Mean error in low uncertainty regions: 0.07
In [493]:
def compute_coverage(y_true, y_pred, scaled_stds):
    """
    Computes the empirical coverage: the fraction of samples where
    the true value is inside the predicted confidence interval.
    
    Parameters:
    - y_true : np.array, ground truth values
    - y_pred : np.array, predicted means
    - scaled_stds : np.array, calibrated standard deviations (after applying k-scaling)

    Returns:
    - coverage : float, percentage of samples inside the confidence interval
    """
    lower_bounds = y_pred - scaled_stds
    upper_bounds = y_pred + scaled_stds

    # Check if the true value lies inside the interval
    inside = (y_true >= lower_bounds) & (y_true <= upper_bounds)

    # Calculate coverage as a percentage
    coverage = np.mean(inside) * 100.0
    return coverage


coverage = compute_coverage(y, y_pred, k_opt * uncertainty)
print(f"Ratio of coverage:{coverage}%")
Ratio of coverage:94.0%

Binning by Mean Predictions

In [494]:
# def calibrate_scaling_adaptive_by_mean(val_pred_means, val_stds, val_targets, num_bins=10, alpha=0.05):
#     """
#     Calibrates uncertainty scaling adaptively based on predicted mean value bins.
#     """
#     # Define bins based on predicted means, not errors.
#     mean_bins = np.linspace(np.min(val_pred_means), np.max(val_pred_means), num_bins + 1)
#     k_values = []

#     for i in range(num_bins):
#         bin_indices = np.where((val_pred_means >= mean_bins[i]) & (val_pred_means < mean_bins[i + 1]))[0]
#         if len(bin_indices) == 0:
#             k_values.append(1.0)  # Default if bin is empty.
#             continue

#         bin_means = val_pred_means[bin_indices]
#         bin_stds = val_stds[bin_indices]
#         bin_targets = val_targets[bin_indices]

#         def coverage_diff(k):
#             lower = bin_means - k * bin_stds
#             upper = bin_means + k * bin_stds
#             coverage = np.mean((bin_targets >= lower) & (bin_targets <= upper))
#             return abs((1 - alpha) - coverage)

#         result = optim.minimize_scalar(coverage_diff, bounds=(0.1, 5.0), method='bounded')
#         k_values.append(result.x)

#     return k_values, mean_bins


# def apply_adaptive_scaling_by_mean(pred_means, pred_stds, k_values, mean_bins):
#     """
#     Applies adaptive scaling based on predicted mean bins.
#     """
#     scaled_stds = np.zeros_like(pred_stds)
#     for i in range(len(k_values)):
#         bin_indices = np.where((pred_means >= mean_bins[i]) & (pred_means < mean_bins[i + 1]))[0]
#         scaled_stds[bin_indices] = pred_stds[bin_indices] * k_values[i]
#     return scaled_stds


# # === Example Usage ===

# # Fit k-values using validation data:
# k_values, mean_bins = calibrate_scaling_adaptive_by_mean(val_pred_means, val_pred_stds, val_targets, num_bins=5, alpha=0.05)

# # Apply learned k-values to a new (or same) dataset:
# scaled_stds = apply_adaptive_scaling_by_mean(val_pred_means, val_pred_stds, k_values, mean_bins)

# print("K values per predicted mean bin:", k_values)
In [495]:
# bin_midpoints = 0.5 * (mean_bins[:-1] + mean_bins[1:])
# plt.plot(bin_midpoints, k_values, marker='o')
# plt.xlabel('Predicted Mean (μ̂)')
# plt.ylabel('Learned Scaling Factor (k)')
# plt.title('Adaptive Scaling Factor per Predicted Mean Bin')
# plt.grid(True, linestyle='--', alpha=0.6)
# plt.show()
In [496]:
# plt.figure(figsize=(8, 6)) # Adjust figure size for publication

# plt.ylim(0, MAX_RUL+1)
# # True Values (more prominent line)
# plt.plot(range(len(y)), y * MAX_RUL, label="True RUL", c="lightseagreen", linewidth=2)

# # Boolean mask: True if point is inside the confidence interval
# inside_interval = (y >= y_pred - scaled_uncertaintity) & (y <= y_pred + scaled_uncertaintity)

# # Assign colors: green if inside, red if outside
# colors_interval = np.where(inside_interval, 'mediumseagreen', 'crimson')

# # Predicted Values (smaller scatter points)
# plt.scatter(range(len(y_pred)), y_pred * MAX_RUL, label="Predicted RUL", c=colors_interval, edgecolor="black", alpha=0.7, zorder=3, s=15) # reduced scatter point size to 15

# # Vertical Lines (improved clarity)
# for i in range(len(y)):
#     plt.plot([i, i], [y_pred[i]*MAX_RUL, y[i]*MAX_RUL], c="gray", linestyle="--", linewidth=0.8, alpha=0.7)

# # Confidence Interval (adjusted opacity, more descriptive label)
# plt.fill_between(
#     range(len(y_pred)),
#     (y_pred - scaled_uncertaintity) * MAX_RUL,
#     np.minimum(y_pred + scaled_uncertaintity, 1) * MAX_RUL,
#     color="salmon",
#     alpha=0.3,  # Reduced opacity
#     label=f"{int(confidence * 100)}% Prediction Interval"  # More accurate label
# )

# # Labels and Title (improved clarity)
# dataset_name = file_name[:-4]
# plt.xlabel("Testing Sample Index") # Removed unnecessary context from label
# plt.ylabel("Remaining Useful Life (RUL)") # Removed unnecessary context from label
# plt.title(f"Comparison of True and Predicted RUL for {dataset_name}") # Removed unnecessary context from title

# # Legend and Grid
# plt.legend()

# # Save and Show
# plt.savefig(f'figures/True_vs_PredictedRUL_simple_calibration_with_bins_{file_name}.png')
# plt.show()
In [497]:
# # Scatter plot with enhancements
# plt.figure(figsize=(8, 6))
# sc = plt.scatter(scaled_uncertaintity, test_errors, c=test_errors, cmap='viridis', alpha=0.7, edgecolor='k')

# # Add colorbar
# cbar = plt.colorbar(sc)
# cbar.set_label('Error Magnitude')

# # Add grid and labels
# plt.grid(alpha=0.3)
# plt.xlabel('Uncertainty', fontsize=12)
# plt.ylabel('Error', fontsize=12)
# plt.title('Uncertainty vs. Error', fontsize=14)

# # Polynomial fit (degree 2 or higher)
# degree = 1  # try 3, 4, etc. for more complex trends
# z = np.polyfit(scaled_uncertaintity, test_errors, degree)
# p = np.poly1d(z)

# # Sort for smooth plotting
# # sorted_idx = np.argsort(scaled_uncertaintity)
# plt.plot(scaled_uncertaintity, p(scaled_uncertaintity),
#           color='red', linestyle='--', label=f'Polynomial Regression (deg {degree})')

# # plt.plot(uncertainty, p(uncertainty), color='red', linestyle='--', label='Trend Line')

# # Add legend
# plt.legend()

# # Show the plot
# plt.tight_layout()
# plt.savefig(f'figures\\UncertaintyVSError_simple_calibration_with_bins_{file_name}.png', dpi=300, bbox_inches='tight')
# plt.show()


# # Compute statistics
# correlation = np.corrcoef(scaled_uncertaintity, test_errors)[0, 1]
# print(f"Correlation between uncertainty and error: {correlation:.2f}")

# # Evaluate performance in high/low uncertainty regions
# threshold = np.percentile(scaled_uncertaintity, 75)  # Upper 25% uncertainty
# high_uncertainty_mask = scaled_uncertaintity > threshold
# low_uncertainty_mask = scaled_uncertaintity <= threshold

# high_uncertainty_error = test_errors[high_uncertainty_mask].mean()
# low_uncertainty_error = test_errors[low_uncertainty_mask].mean()

# print(f"Mean error in high uncertainty regions: {high_uncertainty_error:.2f}")
# print(f"Mean error in low uncertainty regions: {low_uncertainty_error:.2f}")

Callibration with Clustering Algorithm

In [498]:
from scipy import optimize as optim
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
from sklearn.cluster import DBSCAN
In [499]:
X = np.vstack([val_pred_means, val_pred_stds]).T

# Fit Gaussian Mixture Model
num_components = 5
gmm = GaussianMixture(n_components=num_components, random_state=42).fit(X)
cluster_ids = gmm.predict(X)

# Plot clusters
plt.figure(figsize=(8, 6))
scatter = plt.scatter(val_pred_means, val_pred_stds, c=cluster_ids, cmap='tab10', s=40, edgecolor='k', alpha=0.7)

# Plot Gaussian ellipses
def plot_gmm_ellipses(gmm, ax):
    for n in range(gmm.n_components):
        mean = gmm.means_[n]
        cov = gmm.covariances_[n]

        if cov.shape == (2, 2):  # Full covariance matrix
            eigvals, eigvecs = np.linalg.eigh(cov)
            angle = np.degrees(np.arctan2(*eigvecs[:, 1][::-1]))
            width, height = 2 * np.sqrt(eigvals)  # 1-sigma contour
        else:  # Diagonal or spherical covariance
            width, height = 2 * np.sqrt(cov.ravel())
            angle = 0

        ellipse = Ellipse(xy=mean, width=width, height=height, angle=angle,
                          edgecolor='red', facecolor='none', linestyle='--', lw=2)
        ax.add_patch(ellipse)

ax = plt.gca()
plot_gmm_ellipses(gmm, ax)

plt.xlabel("Predicted Mean")
plt.ylabel("Predicted Std")
plt.title("GMM Clustering of Predictions")
plt.legend(*scatter.legend_elements(), title="Cluster")
plt.grid(True)
plt.show()
No description has been provided for this image
In [500]:
def calibrate_scaling_by_clustering(val_pred_means, val_stds, val_targets, num_clusters=3, alpha=0.05):
    """
    Clusters predictions by (mean, std) and calibrates a scaling factor k for each cluster.
    """
    X = np.vstack([val_pred_means, val_stds]).T
    kmeans = KMeans(n_clusters=num_clusters, random_state=42).fit(X)
    cluster_ids = kmeans.predict(X)

    k_values = []

    for cluster in range(num_clusters):
        bin_indices = np.where(cluster_ids == cluster)[0]
        if len(bin_indices) == 0:
            k_values.append(1.0)  # Default scaling if no samples fall here
            continue

        bin_means = val_pred_means[bin_indices]
        bin_stds = val_stds[bin_indices]
        bin_targets = val_targets[bin_indices]

        def coverage_diff(k):
            lower = bin_means - k * bin_stds
            upper = bin_means + k * bin_stds
            coverage = np.mean((bin_targets >= lower) & (bin_targets <= upper))
            return abs((1 - alpha) - coverage)

        result = optim.minimize_scalar(coverage_diff, bounds=(0.1, 5.0), method='bounded')
        k_values.append(result.x)

    return kmeans, k_values

def calibrate_scaling_by_gmm(val_pred_means, val_stds, val_targets, num_components=5, alpha=0.05):
    """
    Calibrates scaling factors k using Gaussian Mixture Model clustering on predicted mean and std.
    """
    X = np.vstack([val_pred_means, val_stds]).T
    gmm = GaussianMixture(n_components=num_components, random_state=42).fit(X)
    cluster_ids = gmm.predict(X)

    k_values = []
    for cluster in range(num_components):
        bin_indices = np.where(cluster_ids == cluster)[0]
        if len(bin_indices) == 0:
            k_values.append(1.0)  # Default scaling for empty clusters
            continue

        bin_means = val_pred_means[bin_indices]
        bin_stds = val_stds[bin_indices]
        bin_targets = val_targets[bin_indices]

        def coverage_diff(k):
            lower = bin_means - k * bin_stds
            upper = bin_means + k * bin_stds
            coverage = np.mean((bin_targets >= lower) & (bin_targets <= upper))
            return abs((1 - alpha) - coverage)

        result = optim.minimize_scalar(coverage_diff, bounds=(0.1, 5.0), method='bounded')
        k_values.append(result.x)

    return gmm, k_values

def calibrate_scaling_by_dbscan(val_pred_means, val_stds, val_targets, eps=0.1, min_samples=10, alpha=0.05):
    """
    Calibrates scaling factors k using DBSCAN clustering on predicted mean and std.
    """
    X = np.vstack([val_pred_means, val_stds]).T
    dbscan = DBSCAN(eps=eps, min_samples=min_samples).fit(X)
    cluster_ids = dbscan.labels_  # -1 means outlier

    unique_clusters = np.unique(cluster_ids)
    k_values = {}

    for cluster in unique_clusters:
        bin_indices = np.where(cluster_ids == cluster)[0]
        if len(bin_indices) == 0:
            continue  # Shouldn't happen

        bin_means = val_pred_means[bin_indices]
        bin_stds = val_stds[bin_indices]
        bin_targets = val_targets[bin_indices]

        def coverage_diff(k):
            lower = bin_means - k * bin_stds
            upper = bin_means + k * bin_stds
            coverage = np.mean((bin_targets >= lower) & (bin_targets <= upper))
            return abs((1 - alpha) - coverage)

        result = optim.minimize_scalar(coverage_diff, bounds=(0.1, 5.0), method='bounded')
        k_values[cluster] = result.x

    return dbscan, k_values


def apply_scaling_by_clustering(pred_means, pred_stds, kmeans, k_values):
    """
    Applies cluster-based scaling using learned k values.
    """
    X = np.vstack([pred_means, pred_stds]).T
    cluster_ids = kmeans.predict(X)

    scaled_stds = np.zeros_like(pred_stds)
    for i, cluster_id in enumerate(cluster_ids):
        scaled_stds[i] = pred_stds[i] * k_values[cluster_id]

    return scaled_stds

def apply_scaling_by_gmm(pred_means, pred_stds, gmm, k_values):
    """
    Applies scaling using GMM-clustered k values.
    """
    X = np.vstack([pred_means, pred_stds]).T
    cluster_ids = gmm.predict(X)

    scaled_stds = np.zeros_like(pred_stds)
    for i, cluster_id in enumerate(cluster_ids):
        scaled_stds[i] = pred_stds[i] * k_values[cluster_id]

    return scaled_stds

def apply_scaling_by_dbscan(pred_means, pred_stds, dbscan, k_values):
    """
    Applies scaling using DBSCAN-clustered k values. Assigns k=1.0 to outliers.
    """
    X = np.vstack([pred_means, pred_stds]).T
    cluster_ids = dbscan.fit_predict(X)  # Predict via re-fit for DBSCAN.

    scaled_stds = np.zeros_like(pred_stds)
    for i, cluster_id in enumerate(cluster_ids):
        k = k_values.get(cluster_id, 1.0)  # Outliers get k=1.0
        scaled_stds[i] = pred_stds[i] * k

    return scaled_stds


# === Example Usage ===

# Fit k-values on validation data:
# kmeans_model, k_values = calibrate_scaling_by_clustering(val_pred_means, val_pred_stds, val_targets, num_clusters=3, alpha=0.05)
gmm_model, k_values = calibrate_scaling_by_gmm(val_pred_means, val_pred_stds, val_targets, num_components=5, alpha=0.05)
# dbscan_model, k_values = calibrate_scaling_by_dbscan(val_pred_means, val_pred_stds, val_targets, eps=0.1, min_samples=10, alpha=0.05)

# Apply to test data:
# scaled_stds = apply_scaling_by_clustering(y_pred, uncertainty, kmeans_model, k_values)
scaled_stds =  apply_scaling_by_gmm(y_pred, uncertainty, gmm_model, k_values)
# scaled_stds =  apply_scaling_by_dbscan(y_pred, uncertainty, dbscan_model, k_values)

print("Cluster-wise learned k values:", k_values)
Cluster-wise learned k values: [0.7750665203530188, 3.128369953115367, 3.8054780887703274, 1.8945177908330837, 3.203146703495102]
In [501]:
plt.figure(figsize=(8, 6)) # Adjust figure size for publication

plt.ylim(0, MAX_RUL+1)
# True Values (more prominent line)
plt.plot(range(len(y)), y * MAX_RUL, label="True RUL", c="lightseagreen", linewidth=2)

# Boolean mask: True if point is inside the confidence interval
inside_interval = (y >= y_pred - scaled_stds) & (y <= y_pred + scaled_stds)

# Assign colors: green if inside, red if outside
colors_interval = np.where(inside_interval, 'mediumseagreen', 'crimson')

# Predicted Values (smaller scatter points)
plt.scatter(range(len(y_pred)), y_pred * MAX_RUL, label="Predicted RUL", c=colors_interval, edgecolor="black", alpha=0.7, zorder=3, s=15) # reduced scatter point size to 15

# Vertical Lines (improved clarity)
for i in range(len(y)):
    plt.plot([i, i], [y_pred[i]*MAX_RUL, y[i]*MAX_RUL], c="gray", linestyle="--", linewidth=0.8, alpha=0.7)

# Confidence Interval (adjusted opacity, more descriptive label)
plt.fill_between(
    range(len(y_pred)),
    (y_pred - 0.95*scaled_stds) * MAX_RUL,
    np.minimum(y_pred + 0.95*scaled_stds, 1) * MAX_RUL,
    color="salmon",
    alpha=0.3,  # Reduced opacity
    label=f"{int(confidence * 100)}% Prediction Interval"  # More accurate label
)

# Labels and Title (improved clarity)
dataset_name = file_name[:-4]
plt.xlabel("Testing Sample Index") # Removed unnecessary context from label
plt.ylabel("Remaining Useful Life (RUL)") # Removed unnecessary context from label
plt.title(f"Comparison of True and Predicted RUL for {dataset_name}") # Removed unnecessary context from title

# Legend and Grid
plt.legend()

# Save and Show
plt.savefig(f'figures/True_vs_PredictedRUL_calibration_with_clustering_bins_{file_name}.png')
plt.show()
No description has been provided for this image
In [502]:
def compute_coverage(y_true, y_pred, scaled_stds):
    """
    Computes the empirical coverage: the fraction of samples where
    the true value is inside the predicted confidence interval.
    
    Parameters:
    - y_true : np.array, ground truth values
    - y_pred : np.array, predicted means
    - scaled_stds : np.array, calibrated standard deviations (after applying k-scaling)

    Returns:
    - coverage : float, percentage of samples inside the confidence interval
    """
    lower_bounds = y_pred - scaled_stds
    upper_bounds = y_pred + scaled_stds

    # Check if the true value lies inside the interval
    inside = (y_true >= lower_bounds) & (y_true <= upper_bounds)

    # Calculate coverage as a percentage
    coverage = np.mean(inside) * 100.0
    return coverage


coverage = compute_coverage(y, y_pred, scaled_stds)
print(f"Ratio of coverage:{coverage}%")
Ratio of coverage:95.0%
In [503]:
# Example usage:
test_errors = np.abs(y - y_pred)

# Scatter plot with enhancements
plt.figure(figsize=(8, 6))
sc = plt.scatter(scaled_stds, test_errors, c=test_errors, cmap='viridis', alpha=0.7, edgecolor='k')

# Add colorbar
cbar = plt.colorbar(sc)
cbar.set_label('Error Magnitude')

# Add grid and labels
plt.grid(alpha=0.3)
plt.xlabel('Uncertainty', fontsize=12)
plt.ylabel('Error', fontsize=12)
plt.title('Uncertainty vs. Error', fontsize=14)

# Polynomial fit (degree 2 or higher)
degree = 1  # try 3, 4, etc. for more complex trends
z = np.polyfit(scaled_stds, test_errors, degree)
p = np.poly1d(z)

# Sort for smooth plotting
# sorted_idx = np.argsort(scaled_uncertaintity)
plt.plot(scaled_stds, p(scaled_stds),
          color='red', linestyle='--', label=f'Polynomial Regression (deg {degree})')

# plt.plot(uncertainty, p(uncertainty), color='red', linestyle='--', label='Trend Line')

# Add legend
plt.legend()

# Show the plot
plt.tight_layout()
plt.savefig(f'figures\\UncertaintyVSError_calibration_with_clustering_bins_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()


# Compute statistics
correlation = np.corrcoef(scaled_stds, test_errors)[0, 1]
print(f"Correlation between uncertainty and error: {correlation:.2f}")

# Evaluate performance in high/low uncertainty regions
threshold = np.percentile(scaled_stds, 75)  # Upper 25% uncertainty
high_uncertainty_mask = scaled_stds > threshold
low_uncertainty_mask = scaled_stds <= threshold

high_uncertainty_error = test_errors[high_uncertainty_mask].mean()
low_uncertainty_error = test_errors[low_uncertainty_mask].mean()

print(f"Mean error in high uncertainty regions: {high_uncertainty_error:.2f}")
print(f"Mean error in low uncertainty regions: {low_uncertainty_error:.2f}")
No description has been provided for this image
Correlation between uncertainty and error: 0.40
Mean error in high uncertainty regions: 0.10
Mean error in low uncertainty regions: 0.05
In [504]:
from scipy.stats import norm

def compute_coverage_at_confidence(y_true, y_pred, scaled_stds, confidence_level):
    """
    Computes the empirical coverage for a given confidence level.
    """
    alpha = 1 - confidence_level
    z = norm.ppf(1 - alpha / 2)  # two-tailed z-score

    lower_bounds = y_pred - z * scaled_stds
    upper_bounds = y_pred + z * scaled_stds

    inside = (y_true >= lower_bounds) & (y_true <= upper_bounds)
    return np.mean(inside)


def compute_ece(conf_levels, emp_coverages):
    """
    Computes Expected Calibration Error (ECE) — L1 distance between confidence and empirical coverage.
    """
    return np.mean(np.abs(np.array(emp_coverages) - np.array(conf_levels)))


def plot_calibration_curve_and_ece(y_true, y_pred, scaled_stds, num_points=20):
    """
    Plots the calibration curve and computes ECE.
    
    y_true: ground truth values
    y_pred: predicted means
    scaled_stds: scaled std deviations (after calibration)
    num_points: number of confidence levels to check between 50% and 99%
    """
    confidence_levels = np.linspace(0.5, 0.99, num_points)
    empirical_coverages = []

    for c in confidence_levels:
        coverage = compute_coverage_at_confidence(y_true, y_pred, scaled_stds, c)
        empirical_coverages.append(coverage)

    # Compute ECE
    ece = compute_ece(confidence_levels, empirical_coverages)

    # Plot
    plt.figure(figsize=(7, 7))
    plt.plot(confidence_levels, empirical_coverages, marker='o', label='Model Calibration')
    plt.plot([0, 1], [0, 1], linestyle='--', color='gray', label='Perfect Calibration')

    plt.xlabel('Expected Confidence Level')
    plt.ylabel('Empirical Coverage')
    plt.title(f'Calibration Curve (ECE = {ece:.4f})')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

    return ece

ece = plot_calibration_curve_and_ece(y, y_pred, scaled_stds)
print(f"Expected Calibration Error (ECE): {ece:.4f}")
No description has been provided for this image
Expected Calibration Error (ECE): 0.1905

Binning by Error Magnitude

In [505]:
# def calibrate_scaling_adaptive(val_pred_means, val_stds, val_targets, num_bins=10, alpha=0.05):
#     """
#     Calibrates uncertainty scaling adaptively based on error magnitude.
#     """
#     errors = np.abs(val_targets - val_pred_means)
#     error_bins = np.linspace(np.min(errors), np.max(errors), num_bins + 1)
#     k_values = []

#     for i in range(num_bins):
#         bin_indices = np.where((errors >= error_bins[i]) & (errors < error_bins[i + 1]))[0]
#         if len(bin_indices) == 0:
#             k_values.append(1.0) #Default value if bin is empty.
#             continue
#         bin_means = val_pred_means[bin_indices]
#         bin_stds = val_stds[bin_indices]
#         bin_targets = val_targets[bin_indices]

#         def coverage_diff(k):
#             lower = bin_means - k * bin_stds
#             upper = bin_means + k * bin_stds
#             coverage = np.mean((bin_targets >= lower) & (bin_targets <= upper))
#             return abs((1 - alpha) - coverage)

#         result = optim.minimize_scalar(coverage_diff, bounds=(0.1, 5.0), method='bounded')
#         k_values.append(result.x)

#     return k_values, error_bins

# def apply_adaptive_scaling(pred_means, pred_stds, errors, k_values, error_bins):
#     """
#     Applies adaptive scaling based on error magnitude.
#     """
#     scaled_stds = np.zeros_like(pred_stds)
#     for i in range(len(k_values)):
#         bin_indices = np.where((errors >= error_bins[i]) & (errors < error_bins[i + 1]))[0]
#         scaled_stds[bin_indices] = pred_stds[bin_indices] * k_values[i]
#     return scaled_stds

# # Example usage:
# k_values, error_bins = calibrate_scaling_adaptive(val_pred_means, val_pred_stds, val_targets, num_bins=10, alpha=0.05)
# val_errors = np.abs(val_targets - val_pred_means)
# scaled_val_stds = apply_adaptive_scaling(val_pred_means, val_pred_stds, errors, k_values, error_bins)

# print("K values per error bin:", k_values)
In [506]:
# # Example usage:
# test_errors = np.abs(y - y_pred)
# scaled_uncertaintity = apply_adaptive_scaling(y, uncertainty, test_errors, k_values, error_bins)
In [507]:
# plt.figure(figsize=(8, 6)) # Adjust figure size for publication

# plt.ylim(0, MAX_RUL+1)
# # True Values (more prominent line)
# plt.plot(range(len(y)), y * MAX_RUL, label="True RUL", c="lightseagreen", linewidth=2)

# # Boolean mask: True if point is inside the confidence interval
# inside_interval = (y >= y_pred - scaled_uncertaintity) & (y <= y_pred + scaled_uncertaintity)

# # Assign colors: green if inside, red if outside
# colors_interval = np.where(inside_interval, 'mediumseagreen', 'crimson')

# # Predicted Values (smaller scatter points)
# plt.scatter(range(len(y_pred)), y_pred * MAX_RUL, label="Predicted RUL", c=colors_interval, edgecolor="black", alpha=0.7, zorder=3, s=15) # reduced scatter point size to 15

# # Vertical Lines (improved clarity)
# for i in range(len(y)):
#     plt.plot([i, i], [y_pred[i]*MAX_RUL, y[i]*MAX_RUL], c="gray", linestyle="--", linewidth=0.8, alpha=0.7)

# # Confidence Interval (adjusted opacity, more descriptive label)
# plt.fill_between(
#     range(len(y_pred)),
#     (y_pred - scaled_uncertaintity) * MAX_RUL,
#     np.minimum(y_pred + scaled_uncertaintity, 1) * MAX_RUL,
#     color="salmon",
#     alpha=0.3,  # Reduced opacity
#     label=f"{int(confidence * 100)}% Prediction Interval"  # More accurate label
# )

# # Labels and Title (improved clarity)
# dataset_name = file_name[:-4]
# plt.xlabel("Testing Sample Index") # Removed unnecessary context from label
# plt.ylabel("Remaining Useful Life (RUL)") # Removed unnecessary context from label
# plt.title(f"Comparison of True and Predicted RUL for {dataset_name}") # Removed unnecessary context from title

# # Legend and Grid
# plt.legend()

# # Save and Show
# plt.savefig(f'figures/True_vs_PredictedRUL_simple_calibration_with_bins_{file_name}.png')
# plt.show()
In [508]:
# # Scatter plot with enhancements
# plt.figure(figsize=(8, 6))
# sc = plt.scatter(scaled_uncertaintity, test_errors, c=test_errors, cmap='viridis', alpha=0.7, edgecolor='k')

# # Add colorbar
# cbar = plt.colorbar(sc)
# cbar.set_label('Error Magnitude')

# # Add grid and labels
# plt.grid(alpha=0.3)
# plt.xlabel('Uncertainty', fontsize=12)
# plt.ylabel('Error', fontsize=12)
# plt.title('Uncertainty vs. Error', fontsize=14)

# # Polynomial fit (degree 2 or higher)
# degree = 1  # try 3, 4, etc. for more complex trends
# z = np.polyfit(scaled_uncertaintity, test_errors, degree)
# p = np.poly1d(z)

# # Sort for smooth plotting
# # sorted_idx = np.argsort(scaled_uncertaintity)
# plt.plot(scaled_uncertaintity, p(scaled_uncertaintity),
#           color='red', linestyle='--', label=f'Polynomial Regression (deg {degree})')

# # plt.plot(uncertainty, p(uncertainty), color='red', linestyle='--', label='Trend Line')

# # Add legend
# plt.legend()

# # Show the plot
# plt.tight_layout()
# plt.savefig(f'figures\\UncertaintyVSError_simple_calibration_with_bins_{file_name}.png', dpi=300, bbox_inches='tight')
# plt.show()


# # Compute statistics
# correlation = np.corrcoef(scaled_uncertaintity, test_errors)[0, 1]
# print(f"Correlation between uncertainty and error: {correlation:.2f}")

# # Evaluate performance in high/low uncertainty regions
# threshold = np.percentile(scaled_uncertaintity, 75)  # Upper 25% uncertainty
# high_uncertainty_mask = scaled_uncertaintity > threshold
# low_uncertainty_mask = scaled_uncertaintity <= threshold

# high_uncertainty_error = test_errors[high_uncertainty_mask].mean()
# low_uncertainty_error = test_errors[low_uncertainty_mask].mean()

# print(f"Mean error in high uncertainty regions: {high_uncertainty_error:.2f}")
# print(f"Mean error in low uncertainty regions: {low_uncertainty_error:.2f}")

Using a Function of Predicted Values to calibrate Uncertaintity

In [509]:
# import scipy.stats as st
# import scipy.optimize as optim

# def calibrate_scaling_by_prediction(val_pred_means, val_stds, val_targets, alpha=0.05, func_type='linear'):
#     """
#     Calibrates uncertainty scaling using a function of the predicted mean (μ).
#     Learns parameters such that the prediction intervals cover (1 - alpha) of the targets.
#     """

#     def compute_coverage(k_values):
#         lower = val_pred_means - k_values * val_stds
#         upper = val_pred_means + k_values * val_stds
#         return np.mean((val_targets >= lower) & (val_targets <= upper))

#     def coverage_diff_linear(params):
#         a, b = params
#         k_values = np.clip(a * val_pred_means + b, 0.01, 10.0)  # prevent negative scaling
#         return abs((1 - alpha) - compute_coverage(k_values))

#     def coverage_diff_quadratic(params):
#         a, b, c = params
#         k_values = np.clip(a * val_pred_means**2 + b * val_pred_means + c, 0.01, 10.0)
#         return abs((1 - alpha) - compute_coverage(k_values))

#     if func_type == 'linear':
#         initial_params = [0.0, 1.0]  # Start from identity scaling
#         result = optim.minimize(
#             coverage_diff_linear,
#             initial_params,
#             method='L-BFGS-B',
#             bounds=[(None, None), (0.1, 5.0)]
#         )
#         optimal_params = result.x
#         print(f"Optimal parameters (linear): {optimal_params}")
#         return lambda pred: np.clip(optimal_params[0] * pred + optimal_params[1], 0.01, 10.0)

#     elif func_type == 'quadratic':
#         initial_params = [0.0, 0.0, 1.0]  # Start from no correction
#         result = optim.minimize(
#             coverage_diff_quadratic,
#             initial_params,
#             method='L-BFGS-B',
#             bounds=[(None, None), (None, None), (0.1, 5.0)]
#         )
#         optimal_params = result.x
#         print(f"Optimal parameters (quadratic): {optimal_params}")
#         return lambda pred: np.clip(optimal_params[0] * pred**2 + optimal_params[1] * pred + optimal_params[2], 0.01, 10.0)

#     else:
#         raise ValueError(f"Unknown func_type: {func_type}")


# def apply_prediction_based_scaling(pred_means, pred_stds, k_function):
#     """
#     Applies the learned scaling function at inference time.
#     """
#     k_values = k_function(pred_means)  # Vectorized, no loop needed
#     scaled_stds = pred_stds * k_values
#     return scaled_stds

# # === EXAMPLE USAGE ===

# # During validation: learn the scaling function
# k_function_pred = calibrate_scaling_by_prediction(
#     val_pred_means, val_stds=val_pred_stds, val_targets=val_targets,
#     alpha=0.05, func_type='quadratic'  #  'linear' or 'quadratic'
# )

# # During inference: apply the learned calibration function
# scaled_test_stds = apply_prediction_based_scaling(y_pred, uncertainty, k_function_pred)

# # Your prediction intervals are now:
# # lower_bound = y_pred - scaled_test_stds
# # upper_bound = y_pred + scaled_test_stds

# print("✅ Calibration function learned and applied.")
In [510]:
# plt.figure(figsize=(8, 6)) # Adjust figure size for publication

# # True Values (more prominent line)
# plt.plot(range(len(y)), y, label="True RUL", c="lightseagreen", linewidth=2)

# # Boolean mask: True if point is inside the confidence interval
# inside_interval = (y >= y_pred - scaled_uncertaintity) & (y <= y_pred + scaled_uncertaintity)

# # Assign colors: green if inside, red if outside
# colors_interval = np.where(inside_interval, 'mediumseagreen', 'crimson')

# # Predicted Values (smaller scatter points)
# plt.scatter(range(len(y_pred)), y_pred, label="Predicted RUL", c=colors_interval, edgecolor="black", alpha=0.7, zorder=3, s=15) # reduced scatter point size to 15

# # Vertical Lines (improved clarity)
# for i in range(len(y)):
#     plt.plot([i, i], [y_pred[i], y[i]], c="gray", linestyle="--", linewidth=0.8, alpha=0.7)

# # Confidence Interval (adjusted opacity, more descriptive label)
# plt.fill_between(
#     range(len(y_pred)),
#     y_pred - scaled_uncertaintity,
#     y_pred + scaled_uncertaintity,
#     color="salmon",
#     alpha=0.3,  # Reduced opacity
#     label=f"{int(confidence * 100)}% Prediction Interval"  # More accurate label
# )

# # Labels and Title (improved clarity)
# dataset_name = file_name[:-4]
# plt.xlabel("Testing Sample Index") # Removed unnecessary context from label
# plt.ylabel("Remaining Useful Life (RUL)") # Removed unnecessary context from label
# plt.title(f"Comparison of True and Predicted RUL for {dataset_name}") # Removed unnecessary context from title

# # Legend and Grid
# plt.legend()

# # Save and Show
# # plt.savefig(f'figures/True_vs_PredictedRUL_{file_name}.png')
# plt.show()
In [511]:
# # Scatter plot with enhancements
# plt.figure(figsize=(8, 6))
# sc = plt.scatter(scaled_uncertaintity, test_errors, c=test_errors, cmap='viridis', alpha=0.7, edgecolor='k')

# # Add colorbar
# cbar = plt.colorbar(sc)
# cbar.set_label('Error Magnitude')

# # Add grid and labels
# plt.grid(alpha=0.3)
# plt.xlabel('Uncertainty', fontsize=12)
# plt.ylabel('Error', fontsize=12)
# plt.title('Uncertainty vs. Error', fontsize=14)

# # Polynomial fit (degree 2 or higher)
# degree = 1  # try 3, 4, etc. for more complex trends
# z = np.polyfit(scaled_uncertaintity, test_errors, degree)
# p = np.poly1d(z)

# # Sort for smooth plotting
# # sorted_idx = np.argsort(scaled_uncertaintity)
# plt.plot(scaled_uncertaintity, p(scaled_uncertaintity),
#           color='red', linestyle='--', label=f'Polynomial Regression (deg {degree})')

# # plt.plot(uncertainty, p(uncertainty), color='red', linestyle='--', label='Trend Line')

# # Add legend
# plt.legend()

# # Show the plot
# plt.tight_layout()
# # plt.savefig(f'figures\\UncertaintyVSError_{file_name}.png', dpi=300, bbox_inches='tight')
# plt.show()


# # Compute statistics
# correlation = np.corrcoef(scaled_uncertaintity, test_errors)[0, 1]
# print(f"Correlation between uncertainty and error: {correlation:.2f}")

# # Evaluate performance in high/low uncertainty regions
# threshold = np.percentile(scaled_uncertaintity, 75)  # Upper 25% uncertainty
# high_uncertainty_mask = scaled_uncertaintity > threshold
# low_uncertainty_mask = scaled_uncertaintity <= threshold

# high_uncertainty_error = test_errors[high_uncertainty_mask].mean()
# low_uncertainty_error = test_errors[low_uncertainty_mask].mean()

# print(f"Mean error in high uncertainty regions: {high_uncertainty_error:.2f}")
# print(f"Mean error in low uncertainty regions: {low_uncertainty_error:.2f}")
In [512]:
# def plot_k_function(k_function, val_pred_means, title="Learned Calibration Function: k(μ)"):
#     """
#     Plots the learned scaling function k(μ) over the range of predicted means.
#     """
#     mu_min = np.min(val_pred_means) - 0.1 * abs(np.min(val_pred_means))
#     mu_max = np.max(val_pred_means) + 0.1 * abs(np.max(val_pred_means))
#     mu_range = np.linspace(mu_min, mu_max, 500)

#     k_values = k_function(mu_range)

#     plt.figure(figsize=(8, 5))
#     plt.plot(mu_range, k_values, label='$k(\\hat{\\mu})$', color='royalblue', linewidth=2)
#     plt.xlabel('Predicted Mean $\\hat{\\mu}$', fontsize=12)
#     plt.ylabel('Scaling Factor $k$', fontsize=12)
#     plt.title(title, fontsize=14)
#     plt.grid(True, linestyle='--', alpha=0.6)
#     plt.legend()
#     plt.show()

# # === Example usage ===

# # Assuming you already ran:
# # k_function_pred = calibrate_scaling_by_prediction(val_pred_means, val_pred_stds, val_targets, alpha=0.05, func_type='quadratic')

# plot_k_function(k_function_pred, val_pred_means)

Intepredability Of the Model¶

SHAP (SHapley Additive exPlanations)

In [513]:
# import shap

# X_list, y_list = [], []

# for batch in valloader:
#     inputs, targets = batch
#     inputs, targets = inputs.to(device).float(), targets.to(device).float()
#     X_list.append(inputs)
#     y_list.append(targets)

# X_np = torch.cat(X_list).cpu().numpy()
# y_np = torch.cat(y_list).cpu().numpy()

# model.eval()

# def model_predict(x_numpy_flat):
#     # Reshape back to original shape before passing to model
#     x_numpy = x_numpy_flat.reshape(-1, 15, 17)
#     x_tensor = torch.from_numpy(x_numpy).float()
#     with torch.no_grad():
#         output = model(x_tensor)
#     return output.view(-1).cpu().numpy()

# # Background subset
# background = X_np[np.random.choice(X_np.shape[0], 100, replace=False)]

# X_flat = X_np.reshape(X_np.shape[0], -1)  # shape: (n_samples, 255)
# background_flat = background.reshape(background.shape[0], -1)  # same here


# # Create explainer
# explainer =  shap.DeepExplainer(model_predict, background_flat)

# # Compute SHAP values for a few samples
# shap_values = explainer.shap_values(X_flat[:10])

# # Plot
# shap.summary_plot(shap_values, X_flat[:10])

Surrogate Model

In [514]:
# from sklearn.tree import DecisionTreeRegressor
# from sklearn.metrics import mean_squared_error

# X_list, y_list, ypred_list =  [], [], []

# for batch in valloader:
#     inputs, targets = batch
#     inputs, targets = inputs.to(device).float(), targets.to(device).float()
#     X_list.append(inputs)
#     y_list.append(targets)
#     ypred_list.append(model(inputs).squeeze())

# X_np = torch.cat(X_list).cpu().numpy()
# y_np = torch.cat(y_list).cpu().numpy()
# ypred_np = torch.cat(ypred_list).detach().cpu().numpy()
# seq, fea = X_np.shape[1], X_np.shape[2]

# X_np = X_np.reshape(X_np.shape[0], -1)

# print(X_np.shape)
# print(y_np.shape)
# print(ypred_np.shape)


# surrogate_model = DecisionTreeRegressor(random_state=42)
# surrogate_model.fit(X_np, ypred_np)
# ypred_surrogate = surrogate_model.predict(X_np)

# mse = mean_squared_error(ypred_surrogate, ypred_np)
# print(f"Mean Squared Error of Surrogate Model: {mse:.4f}")
# feature_importance = surrogate_model.feature_importances_

# # Reshape feature importance back to (time, features)
# importance_matrix = feature_importance.reshape(seq, fea)  # (seq_len, num_features)
# log_importance = np.log10(importance_matrix + 1e-6)

# # Plot heatmap
# plt.figure(figsize=(8, 5))
# sns.heatmap(log_importance.T, cmap="viridis", annot=True, fmt=".2f", cbar=True,
#             xticklabels=[f"T{i+1}" for i in range(30)],
#             yticklabels=[f"F{j+1}" for j in range(4)])

# plt.title("Log-Scaled Feature Importance (Time × Feature)", fontsize=13, weight='bold')
# plt.xlabel("Time Step")
# plt.ylabel("Input Feature")
# plt.tight_layout()
# plt.show()